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Abstract 


One  problem  central  to  the  study  of  slow  deformation  of  the 
ground  is  its  accurate  measurement.  Traditional  techniques, 
repeated  high  precision  surveys,  tilt  meter  observations  and 
strain  meter  observations  suffer  firstly  from  expense, 
difficult  field  procedure  and  low  productivity  and  secondly 
from  the  problem  of  extrapolating  measurements  made  at  a 
single  point  on  the  earth's  surface  to  a  measure  of  average 
behaviour  over  a  larger  area.  It  may  be  possible  to  make 
estimates  of  strain  changes  in  photographic  records  of  the 
ground  at  different  times,  using  commercially  available 
photogrammet r ic  equipment,  combined  with  a  commercially 
available  frequency  analysing  optical  set-up.  Preliminary 
calculations  show  that  it  is  possible  to  detect  a 
differential  ground  displacement  of  1cm,  using  a  simple 
optical  filtering  technique  for  the  analysis  of  the 
photographs . 
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1.  INTRODUCTION 


Prediction  of  occurence  of  large  earthquakes  requires 
careful  monitoring  of  the  accumulation  of  stress  and  strain 
along  major  faults. 

According  to  the  theory  of  plate  tectonics,  the  earth’s 
crust  consists  of  a  number  of  brittle  plates.  Due  to  their 
large  size,  these  plates  will  exhibit  elastic  behaviour,  in 
spite  of  their  relatively  high  shear  strength.  The  plates 
are  driven  across  the  asthenosphere  with  relative  speeds  up 
to  12  cm/year.  Usually,  the  motion  along  a  plate  boundary  is 
continuous,  through  aseismic  creep,  however  portions  of  some 
boundaries  may  become  locked,  with  little  relative  motion 
taking  place.  Stress  accumulates  on  a  locked  fault,  until  it 
reaches  a  critical  value,  known  as  the  failure  stress.  When 
the  failure  stress  is  reached,  the  contact  zone  ruptures, 
causing  an  earthquake,  accompanied  by  a  sudden  displacement, 
which  releases  the  accumulated  stress. 

We  expect  that  due  to  plastic  deformation,  the  fault 
will  be  locked  only  to  a  certain  depth  h,  below  which  the 
motion  will  be  continuous.  Turcotte  and  Spence  (1974)  made 
the  following  simplifications,  which  enabled  them  to  model 
the  strain  accumulation  along  a  locked  strike  slip  fault: 

1.  The  depth  h  is  constant  along  the  fault. 

2.  There  is  zero  friction  below  this  depth. 

3.  The  base  of  each  plate  is  at  a  constant  depth  H,  and  is 
free  to  slide  over  the  asthenosphere. 

Using  these  simplifications,  the  authors  derived  the 
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displacement  field  at  the  surface,  the  mean  stress  on  the 
fault  plane  as  a  function  of  time,  and  the  rate  at  which 
energy  is  stored  in  the  system. 

The  driving  mechanism  of  plates  is  associated  with 
forces  extending  across  the  entire  plate.  Frictional  forces 
acting  on  the  sides  of  the  plates  balance  part  of  the 
driving  forces.  The  accumulation  of  stresses  associated  with 
these  forces  prior  to  an  earthquake  only  extends  a  distance 
q  into  the  plate,  where  q  is  the  typical  dimension  of  the 
fault  break  associated  with  a  strike  slip  earthquake. 

Observations  of  surface  deformations  are  thus  of  vital 
importance  to  the  understanding  of  earthquake  mechanisms. 

Due  to  the  slow  rates  of  strain  accumulation  near  faults, 
measurements  with  a  precision  of  at  least  1  part  per  million 
are  necessary,  if  the  observations  are  to  be  completed  in 
one  or  two  years. 

Satellite  laser-ranging  techniques  have  reached  a  stage 
where  distances  between  two  points  on  the  earth  can  be 
determined  to  within  a  few  centimeters,  and  are  therefore 
exellent  tools  for  determining  motion  on  a  large  scale  (50 
km  or  more).  At  present,  the  main  source  of  error  is  the 
motion  of  the  spacecraft,  which  is  affected  by  gravity 
anomalies . 

For  distances  shorter  than  ten  km,  ground-based 
geodetic  techniques  are  the  principal  means  of  deformation 
monitoring.  In  1973,  Savage  and  Prescott  showed  one  can 
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attain  the  required  accuracy  in  a  single  measurement,  using 
a  Geodolite,  an  electronic  distance  measuring  unit  (EDM). 

A  Geodolite  produces  a  modulated  laser  beam,  which  is 
reflected  back  to  the  instrument  by  a  distant 
ret roref lector .  The  instrument  then  compares  the  modulation 
phases  of  the  incoming  and  outgoing  beams,  which  determines 
the  optical  path  in  terms  of  an  unknown  number  of  whole 
modulation  wavelengths,  plus  a  fraction  of  a  wavelength. 

This  fraction  can  be  precisely  measured,  and  the  number  of 
whole  wavelengths  is  determined  by  repeating  the  measurement 
at  successively  lower  modulation  frequencies.  The  stability 
of  the  modulation  frequency  will  give  the  limit  of  accuracy 
of  the  instrument. 

In  practice  however,  the  accuracy  with  which  we  can 
determine  the  path  length,  or  make  any  optical  measurement, 
will  depend  on  the  accuracy  of  the  measurement  of  the 
refractive  index  of  air,  (n)  through  which  the  beam  passes. 
The  spectrum  of  fluctuations  in  n  shows  significant  power 
from  periods  of  a  few  milliseconds  up  to  very  long  periods, 
measured  in  days  (Slater,  1975).  It  is  reasonable  to  assume 
that  long-period  fluctuations  are  caused  by  weather  changes 
which  affect  very  large  areas,  and  short-period  fluctuations 
are  caused  by  atmospheric  turbulence,  affecting  very  small 
areas  (having  high  spatial  frequency  components).  In  the 
Savage-Prescot t  experiments,  the  high  frequency  fluctuations 
in  n  are  eliminated  by  1-or  10-second  signal  averaging  built 
into  the  geodolite. 


To  eliminate  intermediate  and  long 
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frequency  fluctuations,  it  is  necessary  to  measure  the 
humidity,  temperature  and  pressure  along  the  line  of  sight. 
The  accuracy  with  which  these  measurements  are  made  will 
limit  the  precision  of  the  distance  measurement. 

Savage  and  Prescott  made  these  measurements  by  use  of 
an  aircraft  with  a  probe  mounted  on  each  side  (on  wing 
struts  of  an  airplane  or  skids  of  a  helicopter).  Repeated 
measurements  were  taken  during  the  flight,  to  eliminate 
intermediate  frequency  fluctuations.  Pressure  was  measured 
at  the  two  end  points  of  the  line.  Since  the  accuracy  of  the 
measurement  is  proportional  to  the  stability  of  the 
modulation  frequency,  the  instrument  is  recalibrated  after 
each  line  is  measured.  The  required  agreement  between  the 
modulation  frequency  and  the  frequency  of  a  quartz 
oscillator  was  4 :  10s .  It  was  estimated  that  the  centering  of 
instruments  and  phase  comparison  within  the  Geodolite  cause 
an  error  of  ±1  mm,  and  an  error  in  the  average  pressure, 
temperature  and  vapour  pressure  of  .3  mb,  . 1°C,  and  3  mb 
respectively  cause  an  error  of  .1  ppm  each.  By  staying  in 
continuous  radio  contact  with  the  pilot,  the  flight  path  was 
kept  very  close  to  the  optical  path.  The  airspeed  was  kept 
constant,  and  was  accounted  for  when  determining  the  average 
temperature . 

To  test  the  reproducibility  of  the  data,  the  authors 
resurveyed  a  network  of  30  lines.  The  surveys  were  done 
within  a  three  month  interval,  sc  that  tectonic  movement  did 
not  affect  their  measurements  critically,  but  there  was  a 
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good  meteorological  contrast.  The  differences  in  line 
lengths  were  plotted  against  the  line  length,  and  77%  and 
93%  of  the  points  were  found  to  lie  between  ±c2  and  ±2cr2  , 
where  <x2  ranges  from  4mm  and  10mm  for  lines  1km  and  30km 
respectively . 

The  need  for  an  aircraft  makes  these  measurements 
extremely  expensive,  and  measurements  must  be  done  during 
the  daytime,  as  was  pointed  out  by  Slater  and  Huggett 
(1976).  The  repeatability  of  these  measurements  is  usually 
not  tested  because  of  the  expense.  The  lines  measured  ranged 
from  1  to  35  km  with  standard  deviations  given  by  the 
authors  as  3  mm  to  8  mm  respectively. 

In  1976,  Slater  and  Huggett  developed  and  tested  a 
mult i wavelength  distance  measuring  instrument.  This 
instrument  uses  a  red  He-Ne  laser  (632.9  nm)  and  a  blue 
He~Cd  laser  (441.6  nm)  along  with  a  microwave  source.  The 
difference  in  measured  path  lengths  for  the  two  laser  beams 
gives  an  atmospheric  correction.  Suppose  the  measured  length 
is  L+S,  where  L  is  the  true  length,  and  S  is  the  extra 
length  due  to  the  atmosphere.  S  is  therefore  given  by 


S 


L 

( n- 1 ) d / 


o 


and  the  difference  in  S  obtained  at  two  different 
wavelengths  (red  denoted  by  subscript  one  and  blue  by 
subscript  two)  is  given  by: 
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AS  = 


Sa-S,  = 

y 


(n , -n , ) d / 


A,  ( n , -  1 ) d  7 


(1.1) 


where  A ,  =  ( n 2 -n , ) / ( n , -  1 )  is  only  weakly  dependent  on 
atmospheric  density  and  composition,  so  that  we  can  replace 
it  by  its  average  value  A,  and  take  it  outside  the  integral 
sign  (Bender  and  Owens,  1965).  Therefore,  AS=A1S1.  AS  can 
be  measured,  and  A,  can  be  calculated  from  published  values 
of  the  index  of  refraction  for  the  two  wavelengths,  and  a 
rough  knowledge  of  the  water  vapor  content.  If  AS  is  known 
to  ±.3«,  L  can  be  calculated  to  1  ppm  or  better.  The 
uncertainty  in  water  vapor  pressure  along  the  path  gives  an 
error  of  .1  ppm  per  millibar.  The  index  of  refraction  is 
very  sensitive  to  water  vapor  at  microwave  frequencies, 
therefore  a  third  wavelength  in  the  microwave  range  is  added 
to  the  instrument  to  determine  the  average  vapor  pressure 
along  the  path.  Precision  of  a  few  parts  in  10s,  for  path 
lengths  of  up  to  50  km  and  moderate  temperatures,  should  be 
possible  using  this  instrument.  Field  tests  show  excellent 
results  up  to  distances  of  10  km. 

At  distances  greater  than  10  km,  ambiguity  errors 
arise.  In  this  instrument,  the  modulation  frequency  is 
adjusted  such  that  any  remaining  fraction  of  the  modulation 
wavelength  is  made  zero.  An  error  in  the  measurement  of  k, 
the  number  of  whole  wavelengths,  was  called  ambiguity  error 
by  Slater.  Using  a  single  wavelength,  the  ambiguity  error  is 
resolved  by  slewing  the  modulation  frequency  from  f,  to  f 2 . 
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The  error  in  k  is  given  by  (see  Slater,  1975): 


<fk  =  2k  £f/(  f  2  -f  ,  )  (1.2) 

where  /£  is  the  error  in  £.  The  electronic  stability  of 
their  present  EDM  system  can  determine  the  frequency  to  an 
accuracy  cff /(  f  2 -f  ,  )  =5*  1  0  '  6  .  Therefore  ambiguity  errors  arise 
when  k  becomes  greater  than  10s,  which  corresponds  to  a 
range  of  about  10  km  with  a  modulation  frequency  of  10  cm. 

Another  factor  which  limits  the  range  in  this 
instrument  is  the  signal  to  noise  ratio.  To  allow  successful 
solution  of  the  multi-wavelength  distance  equations,  the 
signal  to  noise  ratio  needs  to  be  high.  This  means  that 
signal  attenuation  due  to  spreading  and  scattering  plays  an 
important  role  in  a  mult i -wavelength  instrument. 

Tectonically  interesting  strains  can  be  observed  using 
small  survey  networks  of  aperture  not  exeeding  2  km 
(Margrave,  1980).  In  fact,  small-aperture  networks  are 
preferred  for  measuring  strain,  to  avoid  averaging  out 
significant  strain  variations  (Savage  and  Prescott,  1973). 

In  his  study  of  the  Peruvian  Andes,  Margrave  used  a  Hewlett 
Packard  3800  EDM.  This  particular  unit  uses  four  different 
modulat i on  f requenc ies  of  an  infrared  beam  to  give  a  reading 
to  within  ±(5  mm  ■+  7  mm/km),  with  a  maximum  range  of  3  km. 
With  a  large  number  of  repetitions  of  distance  measurements, 
they  were  able  to  attain  standard  deviations  comparable  to 
those  obtained  by  Savage  et  al.  in  California.  Atmospheric 
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corrections  were  done  only  at  the  end  points  (at  the 
instrument  and  at  the  reflector),  and  a  linear  relationship 
was  assumed  for  all  points  in  between.  The  instrument  was 
calibrated  before  and  after  the  field  work,  and  the 
modulation  frequency  was  found  to  be  very  stable. 


1 . 1  EDM  Accuracy 

In  the  late  summer  of  1981,  the  HP3800  was  tested  and 
compared  to  a  Sokkisha  Red  1A  EDM.  The  Sokkisha  has  a 
maximum  range  of  only  2  kilometers,  but  is  more  compact, 
easier  to  use  and  can  take  readings  at  a  rate  of 
approximately  once  per  6  seconds.  The  manufacturer's 
specifications  give  the  standard  deviation  as  ±(5  mm  +  5 
mm/km).  The  atmospheric  correction  is  in  steps  of  10  ppm, 
whereas  it  is  continuous  on  the  HP3800.  The  Red  1A  was 
tested  in  the  continuous  read-out  mode,  (figure  1)  at 
approximately  15°C  under  very  windy  and  dry  conditions.  The 
instrument  was  turned  on  at  2:40  pm,  which  corresponds  to 
0.00  seconds  in  the  figure.  Warm-up  time  seems  to  be  about 
100  seconds.  The  instrument  was  accidentally  moved  by  almost 
4  mm  at  1400  seconds,  while  changing  the  battery.  This  would 
not  happen  had  we  used  concrete  pillars  like  those  used  in 
the  microgeodet ic  network  in  Mexico  (Nyland  et  al,  in 
preparation).  From  100  seconds  to  1400  seconds,  the  standard 
deviation  came  to  2.5  mm,  which  is  well  within  the 
manufacturer's  spec i f icat i ons .  and  is  comparable  to 
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.  Performance  of  the  Sokkisha  EDM  under  windy 
0  seconds  corresponds  to  2:40  pm. 
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Margrave's  deviations  which  range  from  1.8  mm  to  4.0  mm 
(Margrave,  1580,  p56).  In  a  comparison  of  the  two  EDM  units 
(figure  2),  the  standard  deviation  was  calculated  using  the 
formula  given  by  Savage  and  Prescott  (1973) 

ct=  (  a  2  +  b2L 2  )-5  (1.3a) 

where  o  is  the  standard  deviation  for  each  instrument.  The 
combined  deviation  of  both  instruments  is  therefore  given  by 
o  2 : 


m2  =  ( 2a  2  +2b 2  L 2  )-s  (1.3  b) 

For  an  ideal  normal  distribution,  68%  and  95%  of  the  points 
lie  between  ±c  and  ±2a  respectively,  therefore  we  can  take 
a=3  mm  and  b=3*10~6.  These  values  give  70%  and  97%  of  the 
points  between  ±o2  and  ±2a2  . 

There  seems  to  be  a  slight  bias  toward  the  positive, 
which  means  that  the  Sokkisha  gave  cons i stant l.y  higher 
values  for  the  same  distance  than  the  HP3800.  This  was 
probably  due  to  the  instrument  constant,  which  was  later 
found  to  need  slight  adjustment.  To  test  this,  two 
reflectors  were  placed  about  100  meters  apart,  with  the 
instrument  roughly  half  way  between  them  and  the  distance  to 
each  reflector  was  measured  and  added  to  give  the  total 
distance.  One  of  the  reflectors  was  then  replaced  by  the 
instrument,  and  the  total  distance  was  measured  directly. 


1 1 


L  (meters) 


Figure  2....  Difference  between  distance  measurements  of 
the  SOKKISHA  and  HP3800  infrared  EDM  units.  70%  of  the 
points  lie  between  ±o"2  ;  97%  lie  between  ±2 c"2 


1  2 


The  difference  between  these  two  values  is  called  the 
instrument  constant. 

Figures  3  and  4  show  the  percentage  of  return  signal  as 
a  function  of  area  of  reflector  surface  (figure  3),  and  as  a 
function  of  distance  from  the  center  of  the  reflector 
(figure  4).  The  maximum  return  signal  is  dependent  on  the 
distance  of  the  instrument  from  the  reflector,  and  on  the 
weather  conditions.  From  figure  4,  the  spreading  of  the  beam 
can  be  calculated  to  be  about  10" 3  radians.  Due  to  the 
spreading  of  the  infrared  beam  and  the  use  of  corner  cube 
reflector,  the  pointing  of  the  instrument  at  the  reflector 
is  not  critical;  a  good  reading  can  be  taken  with  a  return 
signal  of  25%  or  more.  The  background  level  was  about  15%  in 
both  cases  (fig.  3  and  4);  readings  were  taken  under  hot  and 
windy  weather  conditions.  The  three  corner  cube  reflectors 
(5  cm  diameter  each)  were  set  up  in  a  triangular 
configuration  with  the  centers  of  the  prisms  separated  by 
6cm .  The  effective  radius  of  the  reflector  was  about  5  cm. 

Even  with  stable  monumen tat i on , repea ted  surveys  of 
geodetic  networks  are  expensive  and  require  a  high 
commitment  of  manpower.  Under  normal  c i rcums tances ,  one 
should  allow  2-3  weeks  for  a  survey  of  a  typical  network 
(35-40  distances  and  approx.  80  angles)  (Margrave , 1 980 ) . 

I  propose  an  alternative  observation  method.  If  the 
position  of  discrete  points  on  the  surface  of  the  earth  is 
recorded  photographically  at  different  times,  a  comparison 
of  the  two  records  by  an  optical  Fourier  technique  might 


Figure  3....  Fall-off  of  return  signal  with  decreasing 
reflector  surface  area  for  three  different  distances  for  tn 
HP3800 
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cd/rd 


Figure  4....  The  return  signal  reaching  the  instrument 
pointed  a  distance  D  away  from  the  reflector  of  radius  Rf 
1290  meters  away.  Spreading  of  the  beam  is  therefore  about 
.00  1  radians  . 
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serve  to  determine  the  strain  which  has  accumulated  during 
that  time. 

One  point  on  the  ground  will  be  taken  as  the  reference 
point.  Limitations  on  the  precision  of  camera  alignment  will 
result  in  slight  displacement  of  and  rotation  about  this 
point.  The  method  is  thus  capable  of  measuring  accumulated 
strain,  but  not  the  actual  displacements.  The  attainable 
accuracy  is  limited  by  the  atmosphere,  as  in  the  previously 
discussed  methods.  Atmospheric  conditions  for  desert  areas, 
and  measurements  of  atmospheric  turbulence  are  discussed  in 
detail  by  Walters  and  Kunkel  (1981). 

It  is  the  aim  of  this  work  to  find  the  theoretical 
feasibility  of  this  method. 


2.  IMAGING  SYSTEMS 


To  understand  the  approach  we  will  take  in  analyzing 
our  two  photographs,  we  need  to  first  look  at  imaging 
systems  and  how  they  work.  In  particular,  we  want  to  look  at 
the  part  of  imaging  systems  theory  called  frequency  analysis 
of  optical  systems,  or  Fourier  optics.  In  this,  and  all 
following  chapters,  'frequency'  is  taken  to  mean  the  spatial 
frequency.  When  speaking  about  the  light  source,  I  will 
refer  to  its  wavelength  only,  to  avoid  confusion.  It  is 
beyond  the  scope  of  this  chapter  to  give  more  than  a 
superficial  treatment  of  Fourier  optics;  a  number  of  good 
books  are  available  on  the  subject,  (eg,  Goodman,  1968) 

As  an  example,  consider  a  simple  imaging  system, 
consisting  of  a  point  source,  a  camera  and  film.  In  the 
resultant  photograph,  the  point  will  be  spread  by  an  amount 
determined  by  the  resolution  of  the  film  and  the  size  of  the 
camera  aperture.  That  is,  both  the  camera  and  the  film  have 
a  point  spread  function  associated  with  them. 

In  any  real  camera,  lens  aberrations,  and  errors  due  to 
the  camera  frame  cannot  be  made  exactly  zero.  These 
aberrations  can  be  minimized  in  a  compound  lens  consisting 
entirely  of  spherical  surface  lenses,  by  use  of  elaborate 
computer  programs  which  find  optimum  configurations,  and 
appropriate  multi-layer  coatings.  In  good  terrestrial 
cameras,  the  aberration  errors  are  small,  and  some  of  these 
errors  will  cancel  out  when  we  compare  two  photographs  taken 
with  the  same  camera  and  lens,  since  we  are  only  interested 
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in  relative  displacements. 


2 . 1  Lenses 

The  most  important  components  of  optical  systems  are 
lenses.  Consider  a  thin  converging  spherical  lens  (a  lens  is 
said  to  be  thin  if  its  radii  of  curvature  are  much  greater 
than  its  thickness).  A  lens  is  made  of  material  optically 
denser  than  air,  and  hence  it  will  delay  a  light  ray  by  an 
amount  proportional  to  the  thickness  of  the  lens  at  that 
point.  Assume  we  have  a  spherical  lens  of  thickness  T(x,y), 
with  a  maximum  thickness  T0.  If  g(x,y)  is  the  complex 
amplitude  distribution  immediately  in  front  of  the  lens,  and 
g'(x,v)  is  the  distribution  immediately  behind  the  lens, 
(figure  5)  then 


g '  =  t  g  (2.1) 

where  t  is  the  phase  transformation  and  is  given  by: 

t=exp(ikT0)exp(ik(n-1)T(x,y))  (2.2) 

In  this  equation  the  term  containing  nT  is  due  to  the  phase 
delay  in  the  lens,  and  the  term  containing  (T0-T)  is  due  to 
the  delay  by  the  space  adjacent  to  the  lens,  k  is  the  wave 
number,  given  by  k  =  2n/\ 
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Figure  5  .  .  .  .  Object 
g,is  placed  in  front 
focal  length  f.  The 


with  complex  amplitude  distribution 
of  a  convex  thin  spherical  lens  of 
resulting  amplitude  distribution  is  g 
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From  the  geometry  of  a  lens  with  spherical  surfaces 
with  radii  of  curvature  R,  and  R2  we  can  calculate  the 
thickness  of  the  lens  at  any  point  (x, y): 

T=T  o - ( x  2  +y  2 ) /2  (1/R,  -  1/R2) 

=  T  o  -  ( x  2  +v 2 ) /( 2f ( n- 1 ) )  (2.3) 

where  we  used  the  lensmaker's  formula  in  the  last  step.  The 
phase  transformation  will  therefore  be  (if  we  substitute 
(2.3)  into  (2.2)  and  drop  the  constant  phase  delay, 
exp[ i knT  0 ] ) : 


t  =  exp[-ik(x2+y2 ) / 2 f  3 


(2.4) 


The  complex  amplitude  distribution  behind  the  lens  is 
therefore  given  by: 

g ' ( x , y ) =g ( x , y ) exp[ - i k ( x  2  +y 2 ) /2f ]  (2.5) 

One  of  the  most  useful  properties  of  a  converging  lens 
is  its  ability  to  perform  a  two  dimensional  Fourier 
transform.  The  back  focal  plane  of  a  converging  lens 
contains  the  Fourier  transform  of  a  coherently  illuminated 
object  placed  in  the  front  focal  plane  (proof  is  given  in 
the  Appendix).  The  back  focal  plane  is  therefore  called  the 
transform  plane,  or  frequency  plane,  and  the  front  focal 
plane  is  usually  referred  to  as  the  input  plane. 
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In  interferometry,  we  often  talk  about  phase 
differences  of  a  given  number  of  wavelengths.  Since  fully 
monochromatic  light  is  impossible  to  attain,  we  cannot  have 
a  path  difference  of  one  wavelength  precisely  for  all  the 
wavelengths  in  the  beam.  The  degree  of  monochromaticity  of 
the  light  is  called  temporal  coherence. 

Because  mathematical  point  sources  are  excludable  in 
principle,  rays  reach  a  point  in  an  interference  pattern 
from  all  points  in  the  source  "point".  Perfect  matching  is 
therefore  impossible.  The  degree  of  approximation  of  a 
source  to  a  point  source  is  called  spatial  coherence. 

Before  the  invention  of  lasers,  good  coherence  was 
achieved  by  spectral  filtering  and  by  making  the  source  very 
small  (resulting  in  considerable  light  loss).  Lasers  can 
only  amplify  photons  if  their  wavelengths  and  directions  of 
travel  lie  in  certain  narrow  bands,  and  therefore  good 
coherence  is  achieved  without  light  loss. 

The  effects  of  coherence  on  imaging  systems  have  been 
discussed  by  P.  Considine  (1966),  and  will  be  summarized  in 
the  next  section. 


2.2  Effects  of  coherence  on  imaging  systems 

The  highest  frequency  passed  by  a  lens  depends  on  its 
aperture.  Suppose  we  are  trying  to  image  a  sharp  edge,  (such 
as  that  of  a  razor  blade)  containing  very  high  frequency 
components.  A  sketch  of  intensity  versus  distance 


perpendicular  to  the  edge  was  made  (figure  6).  An 
examination  of  a  photograph  of  a  coherently  or  incoherently 
illuminated  edge  will  show  this  type  of  an  intensity 
distribution.  The  same  distribution  can  be  derived 
theoretically  (see  Considine,  1966).  From  the  sketch,  two 
things  are  immediately  apparent: 

1.  The  coherent  illumination  produces  ringing  or  intensity 
oscillations  about  I0  (the  Fresnel  diffraction  of  an 
edge ) . 

2.  On  the  basis  that  the  edge  is  located  at  the 
half-intensity  point,  the  coherently-illuminated  edge  is 
shifted,  (see  Considine,  1966) 

The  two  effects  are  a  result  of  the  difference  between 
transfer  functions  of  coherent  and  incoherent  imaging 
systems.  That  is,  the  transfer  funcion  of  an  optical  system 
depends  on  the  type  of  illumination.  In  the  frequency  plane, 
a  coherently  illuminated  lens  has  a  sharply  cut-off  transfer 
function,  whereas  the  transfer  function  of  an  incoherently 
illuminated  lens  decreases  gradually  with  frequency.  If  the 
cut-off  frequency  is  lower  than  the  highest  frequency  of  the 
edge,  ringing  will  occur  in  the  coherent  system.  This  is 
important,  since  ringing  greatly  influences  the  resolution. 
Considine  shows  that  a  lens  which  can  resolve  114  lines/mm 
when  incoherently  illuminated,  only  resolves  36  lines/mm 
when  a  coherent  light  source  is  used,  because  the  ringing  of 
one  line  starts  to  overlap  the  next  closest  line. 


t-b  *~rj 
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igure  6....  Image  of  a  sharp  edge:  Intensity  as  a  function 
distance  perpendicular  to  the  edge.  Dashed  line 
represents  coherent  illumination;  solid  line  represents 
incoherent  illumination  of  the  edge 
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Another  effect  of  coherent  illumination  is  speckling, 
which  is  not  a  property  of  the  light  beam  itself,  but  is 
present  only  when  phase  variations  are  introduced  into  the 
wavefront  by  transmission  through  an  optically  rough  surface 
such  as  a  scotch  tape,  or  by  reflection  from  a  rough  surface 
such  as  a  painted  wall.  The  size  of  the  speckles  depends  on 
the  wave  number,  and  they  can  thus  be  made  small  enough  to 
be  useful  in  certain  types  of  experiments  (Khetan  and 
Chiang,1376)  which  will  be  discussed  later. 

There  is  no  sign  of  edge  ringing  within  the  speckle 
pattern.  A  simple  experiment  can  show  that  the  absence  of 
edge  ringing  is  not  caused  by  loss  of  coherence: 

Passing  a  laser  beam  through  a  diffuser  such  as  a 
scotch  tape  produces  a  speckle  pattern.  Placing  a  double 
slit  in  contact  with  the  diffuser  produces  clear 
interference  fringes  across  the  speckle  pattern.  The  field 
thus  remained  coherent  after  having  been  diffused. 


2.3  Optical  Processing 

If  g,  is  the  input  function  and  g2  is  the  corresponding 
output  for  a  linear  system, 


g  2  (  x  '  ,  y  '  )  = 


no 


g 


( x , y )p( x ’ , y ' ;x,y) dxdy 


(2.6) 


-00 


where  p  is  the  impulse  response,  or  the  point  spread 
function.  If  the  system  is  space- invar iant , 


then  p  depends 
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onlv  on  the  distances  (x'-x)  and  (v 


g :  ( x  ’  f 


y  ’ )  =  ^  g ,  ( x ,  y )  p  ( x 


-X  ,y 


-y ) dxay 


(2.7) 


which  is  iust  the  two  dimensional  convolution: 


(2.8) 


Taxing  the  Fourier  transform,  we  get: 


G  2 ( u , v ) =G ,  ( u , v )  ? ( u , v  ) 


(2.9) 


where  ?  is  the  Fourier  transform  of  p,  and  is  called  the 
transfer  function  of  the  system. 

Consider  a  plane  wave  illuminating  an  input  placed  at 
plane  P, .  (figure  7)  A  thin  spnerical  lens  L2  is  placed  a 
distance  f,  from  where  f ,  is  the  focal  length  of  the 

lens.  In  the  back  focal  plane  of  the  lens  will  appear  the 
Fourier  transform  of  che  object  .  This  plane  is  called  the 
transform  plane,  and  is  denoted  by  Pa .  We  can  Fourier 
transform  this  plane  by  the  use  of  another  lens.  The  output 
will  appear  in  plane  P3 ,  the  image  plane.  The  image  is 
inverted,  because  we  are  performing  two  consecutive 
transforms,  rather  than  a  Fourier  transform  followed  by  an 
inverse  transform.  We  therefore  choose  the  coordinates  in  P3 
in  reversed  directions,  as  denoted  by  the  negacive  sign  in 
figure  7.  The  magnification  is  given  by  the  ratio  of  the 


Figure  7. 


•  •  • 


Arrangement  of  an  optical  Fourier  analyser 
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focal  lengths  of  the  two  lenses,  and  is  usually  chosen  to  be 
unity.  The  plane  P2  contains  the  Fourier  transform  of  the 
input,  with  the  spatial  frequencies  related  to  the  spatial 
coordinates  by: 


x  2  =  f  A  u 


(2.10) 


y 2=f Av 

By  placing  appropriate  filters  in  the  transform  plane,  we 
can  therefore  alter  the  image  in  the  output  plane.  For 
example,  to  block  out  high  frequency  components  of  the 
input,  we  would  place  a  small  aperture  at  the  transform 
plane.  The  highest  frequency  passed  by  a  given  aperture  is 
given  by  equations  2.10. 

The  first  attempts  to  manipulate  the  Fourier  spectrum 
of  an  image  were  done  in  1893  by  Abbe,  and  a  few  years  later 
by  Porter.  The  Abbe-Porter  experiments  are  a  good 
demonstrati  on  of  the  principle  of  Fourier  analysis.  A  fine 
wire  mesh  is  placed  in  a  collimated  coherent  beam,  in  the 
front  focal  plane  of  a  lens.  The  back  focal  plane  will 
contain  a  rectangular  array  of  dots,  the  Fourier  transform 
of  the  mesh.  With  a  narrow  horizontal  slit,  passing  only  one 
row  of  spectral  components,  the  image  will  contain  only  the 
vertical  components  of  the  wire  mesh  input. 

Most  useful  filters  are  harder  to  fabricate,  and  until 
the  development  of  digital-optical  converters,  they  could 
only  be  approximated.  These  filters  are  used  to  attenuate 


undesired  effects  such  as  linear  smear,  blur  or  atmospheric 
effects  in  photographs.  It  is  hoped  that  the  transfer 
function  of  the  compensating  filter  will  remove  the  transfer 
function  of  the  undesired  effect,  therefore  the  filter  is 
given  by  (if  we  ignore  noise): 


F ( y )  =  1  /  p(?) 


(2.11) 


where  P(y)  is  the  Fourier  transform  of  the  point  spread 


function  of  the  effect  we  want  to  remove.  For  example,  the 
point  spread  for  linear  smear  is  a  rectangle  function,  and 
the  transfer  function  is  simply  its  Fourier  transform,  which 
is  the  Sine  function. 

If  additive  noise  is  present,  however,  the  filter  given 
by  equation  2.11  could  degrade  the  photograph  instead  of 
improving  it,  depending  on  the  signal  to  noise  ratio.  To 
include  the  noise,  we  start  with  the  equation  for  a  linear 
system  (eqn.  2.6).  In  this  example,  our  system  consists  of 
an  object,  atmosphere,  a  camera  and  film.  The  output 
o(x',y’)  is  related  to  the  input  i(x,y)  by 


(2.12  a  ) 


0 ( u , v ) =1 (u,v)P(u,v) 


(2.12  b) 


P  contains  all  lens  aberrations  and  conditions  which  may 
degrade  the  image.  The  overall  transfer  function  is 
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therefore  given  by: 


P  =P,PLPnPr 


A  rL  rn 


(2.13) 


where : 

P„  =atmospheric  transfer  function, 

PL  =  transfer  function  of  the  lens 

Pn  =  transfer  function  of  motion  during  exposure 

PF  =  transfer  function  due  to  the  film. 

We  now  complete  the  model  by  considering  noise.  The 
noise,  which  we  will  assume  to  be  additive,  consists  of  all 
the  random  processes  present,  such  as  dirt,  dust, 
discontinuities,  film  grain  and  laser  speckle.  The  observed 
image  is  therefore 


(2.14) 


o(x,y)=i(x,y)*p(x,y)+n(x,y) 


If  we  take  i(x,y)  to  be  the  estimate  of  i(x,y),  then  the 


error  function  e(x,y)  is  just  i(x,y)  -  i(x,y),  and  has  the 


power  spectrum  |E(u,v)|2.  The  mean  square  error  will  be 
given  by: 


dudv 


where  W  is  a  frequency  weighting  function.  For  example,  for 
frequencies  which  we  consider  unimportant,  we  can  put  W 
equal  to  zero. 
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If  we  assume  that  the  noise  is  not  correlated  with  the 
signal,  then  for  a  linear  system,  minimizing  the  mean  square 
error  gives  (see  Considine  and  Gonsalves,  1978): 


F(/J  = 


1 

0(37) 

N(y) 

P(p) 

067) 

N(7) 

2 

+ 

P(y; 

-2 

(2.15) 


which  is  knov7n  as  the  Wiener  filter.  It  was  first  designed 
to  improve  atmosphere-degraded  images.  The  estimate  of  the 
input,  'l  is  0  ( u  ,  v  )  »F  (  u  ,  v )  .  It  can  be  seen  that  if  the  signal 
to  noise  ratio  is  very  large,  this  reduces  to  the  inverse 
filter,  1/P.  Thus  for  a  system  with  infinite  signal  to  noise 
ratio,  we  get: 

1  =  OF  =  I  P/P  =  I  (2.16) 


2.4  Atmosphere 

In  this  section,  we  assume  that  in  our  system,  there  is 
no  limitation  on  size  or  quality  of  the  lens,  or  resolution 
of  the  film.  The  light  beam  we  are  recording  will  be 
distorted  by  inhomogeneities  in  the  refractive  index  of  the 
atmosphere  through  which  it  travels.  D.L.  Fried  (1966) 
calculated  the  effects  of  the  atmosphere  on  a  point  source 
as  a  function  of  altitude,  using  a  model  of  vertical 
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distribution  of  the  index  of  refraction  structure  parameter, 
C*  based  on  experimental  dara.  He  found  that  at  an  altitude 
h , 


C2  =4.2*10“  1  4  h“,/3  exp(-h/h0) 


(2.17) 


where  h0  was  found  to  be  3200  m.  C*  is  a  coefficient  which 
measures  the  strength  of  the  refractive  index  turbulence,  a 
function  of  both  time  and  space.  It  is  usual  to  measure  the 
temperature  structure  parameter  Cj  ,  which  is  directly 
proportional  to  C2M  .  Another  parameter  which  is  often  used  to 
give  a  measure  of  turbulence  is  the  transverse  coherence 
length,  r0.  The  relation  between  C*  and  r0  is: 


c 


r 


a 


[  C^(z)dz]“3/5 
J 
0 


(2.18) 


For  a  long  exposure,  the  maximum  resolution  R.  has  been 
defined  as  the  limiting  value  of  the  resolution  R,  as  the 
diameter  of  the  lens  becomes  infinitely  large.  In  terms  of 
the  altitude  z,  Fried  gets: 


R.  =  ( n/4 ) ( r 0 /\  z ) 2  (cycles/meter)2  (2.19) 


At  an  altitude  z,  there  is  a  minimum  size  of  an  object  which 
we  can  see.  The  maximum  spatial  frequency  is  given  by  the 
square  root  of  the  resolution,  and  is  in  units  of  line  pairs 
per  meter.  The  minimum  resolvable  length  is  : 
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/  =  1  /  2  (  R ) ,s  (2.20) 

R  approaches  its  limiting  value  asymptotically  as  the  lens 
diameter  D  becomes  very  large.  When  D=2r0,  R  is  almost  70% 
of  R+ .  When  D  equals  r0,  R  is  approximately  .45R+. 

The  transverse  coherence  length  r0  has  been  measured 
accurately  for  mountain  and  desert  areas  by  Walters  &  Kunkel 
(1981).  During  the  night,  r0  in  the  desert  is  fairly 
constant  until  sunrise,  with  a  mean  value  of  50mm.  Optical 
turbulence  is  lowest  (r0  is  at  a  maximum)  when  the 
temperature  of  the  ground  is  equal  to  the  temperature  of  the 
air.  This  occurs  just  after  sunrise,  and  just  before  sunset. 


2.5  Optical  monitoring  of  surface  deformations 

The  method  proposed  here  is  a  variation  of  the  idea  of 
Vacquier  and  Whiteman  (1973),  who  used  a  camera  of  10-meter 
focal  length  to  photograph  two  pairs  of  lights,  one  on  each 
side  of  an  active  fault.  Average  displacement  over  the  past 
5  million  years  on  this  portion  of  the  fault  is  6  cm/vear. 
The  separations  of  the  dots  on  the  photographic  plates  were 
measured  using  a  10-power  microscope  with  a  micrometer 
screw.  If  there  is  a  horizontal  temperature  gradient,  a 
systematic  error  would  be  introduced.  This  error  would  be 
canceled  out  by  another  camera,  pointed  in  the  opposite 
direction.  This  second  station  was  never  built,  mainly 
because  of  the  expense,  and  because  the  authors  felt  that  it 
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would  not  significantly  improve  the  precision. 

Over  a  two  year  period,  the  expected  12  cm  displacement 
was  not  observed.  The  authors  believe  this  is  because  the 
fault  is  locked,  and  over  the  distance  of  25  km  (the 
distance  between  the  camera  and  the  far  pair  of  lights)  the 
strain  is  linear.  Since  they  were  measuring  parallax,  only  a 
non-linear  strain  could  be  detected. 

Even  though  their  experiment  did  not  give  the  expected 
results,  much  can  still  be  learned  from  their  data  (p  362, 
JGR,  vol.78,  no. 5).  There  is  a  great  variation  in  the 
quality  of  the  photographs.  The  size  of  the  spots  can  vary 
from  .15  to  .5  mm.  The  photographs  were  taken  in  twenty 
minute  intervals,  in  which  the  size  of  the  spots  could 
change  by  a  factor  of  2  or  more. 

Two  pairs  of  dots  appear  on  each  photograph:  the  inner 
pair  being  the  image  of  the  far  lights  (25  km),  and  the 
outer  pair  the  near  lights  (8  km).  The  inner  pair  may  appear 
above  or  below  the  outer  pair,  depending  on  the  vertical 
temperature  gradient  during  the  exposure.  The  vertical 
displacement  of  the  inner  dots  varies  from  about  .2  mm  to 
-.1  mm.  Horizontal  temperature  gradient  cannot  be  observed 
from  these  photographs. 

The  expected  shift  of  the  inner  pair  of  spots  is  .05  mm 
over  the  two  year  period.  Averaging  over  a  large  number  of 
photographs,  the  authors  claim  to  be  able  to  detect  a  shift 
of  .002  mm  (corresponding  to  a  ground  movement  along  the 


fault  of  4  mm). 
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Our  method  will  differ  in  two  ways.  One,  we  will  use  a 
random  array  of  light  sources  or  reflectors,  and  two,  we 
will  process  the  data  using  a  double  exposure 
interferometric  technique.  The  area  of  study  (approximately 
1  square  kilometer)  will  be  an  active  fault  in  a  desert  area 
in  northern  Mexico  (Darby  et  al.,  1981).  The  motion  along 
the  fault  is  almost  entirely  horizontal,  and  expected 
maximum  yearly  displacements  are  1  to  2  cm.  An  array  of 
precisely  positioned  reflectors  in  the  area  will  be  imaged 
on  a  10  by  10  cm  high  resolution  film,  on  an  Estar  base  or 
glass  plate  for  good  stability.  Holographic  films  with  a 
resolution  of  500  to  2000  lines/mm  are  commercially 
available  on  long  lasting  stable  bases.  For  more  discussion 
on  recording  media,  see  Goodman  (1968). 

The  recording  of  the  lights  is  thus  the  same  as  the 
Vacquier  and  Whiteman  experiment  on  a  smaller  scale.  Where 
their  line  of  sight  was  25  km,  ours  will  be  about  1  km. 

Their  10  meter  focal  length  lens  has  a  resolution  of 
approximately  25  microns,  and  the  displacement  they  were 
looking  for  was  50  microns.  Our  lens  will  have  a  focal 
length  5  to  10  cm,  with  a  resolution  .5  to  1  micron,  and  the 
displacement  we  are  looking  for  is  1  to  2  microns.  Our 
system  will  be  manually  operated  whereas  theirs  was 
automatic,  but  Vacquier  and  Whiteman  ran  into  some  problems 
with  their  automatic  system,  and  most  of  their  useful  data 
was  taken  manually.  Photographs  should  be  taken 
approximately  every  twenty  minutes,  because  atmospheric 
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conditions  can  change  significantly  in  this  time  interval, 
and  the  size  of  the  images  of  the  reflectors  will  vary.  When 
we  compare  two  transparencies  separated  in  time  by  one  or 
two  years,  we  will  choose  the  two  in  which  the  size  of  the 
reflector  images  is  similar. 

It  is  assumed  that  an  area  10  meters  in  length  will 
have  a  uniform  displacement  rate,  and  that  it  is  sufficient 
to  use  an  array  of  reflectors  with  a  separation  of  about  10 
meters.  The  transparencies  will  be  compared  using  a  double 
exposure  interferometry  technique.  The  attainable  accuracy 
will  be  discussed  in  chapter  4. 

The  position  of  the  reflectors  is  not  critical,  but  a 
line  perpendicular  to  the  fault  and  passing  through  a  point 
directly  beneath  the  camera  (referred  to  as  the  main  axis 
from  now  on)  will  have  the  least  error,  since  the  vertical 
index  of  refraction  gradient  only  affects  precision  on  the 
component  perpendicular  to  the  fault.  The  reflectors  should 
therefore  be  concentrated  near  the  main  axis,  as  shown  in 
the  lower  part  of  figure  8.  (Black  squares  represent 
reflectors).  When  in  the  field,  it  is  usually  impossible  to 
make  a  straight  row  of  lights  which  does  not  deviate  from 
the  main  axis.  It  is  therefore  useful  to  estimate  the  error 
introduced  by  the  vertical  refraction  gradient.  A  camera  on 
top  of  a  mountain  at  a  height  h,  is  pointed  at  point  P  below 
(where  P  is  the  position  of  the  reflector  which  deviates 
most  from  the  main  axis).  The  light  rays  will  follow  a  bent 
path,  and  the  image  of  P  will  seem  to  originate  at  P' 
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( f i gure  8 )  . 

The  ground  in  a  desert  area  at  night  cools  off  very 
rapidly.  Dense  air  near  the  ground  is  cooled,  and  thus  tries 
to  maintain  its  position  near  the  ground.  Thus  atmospheric 
turbulence  is  near  minimum  during  the  desert  night.  A 
detailed  discussion  on  measurements  of  atmospheric 
turbulence  in  a  desert  area  can  be  found  in  a  paper  by 
Walters  and  Kunkel  (1981). 

The  assumption  I  make  is  that  the  index  of  refraction 
varies  linearly  with  altitude  z.  That  is,  n(z)=c+kz.  From 
Snell's  law  we  have  that  n(z)sin@(z)  is  a  constant.  The 
distance  s  is  thus  given  by: 


s=  dz  tan6 ( z ) 


0 


=(a/k) [cosh" 1 ( (c+kh)/a)-cosh‘ ’ (c/a) ] 


(2.21  ) 


where  a=sin/  ,  and  9  is  the  angle  at  any  point  along  the 
curved  path  to  the  vertical.  As  a  specific  example,  suppose 
the  vertical  temperature  gradient  is  5  degrees  per  100 
meters  of  altitude.  This  value  is  extremely  high;  in  most 
cases,  the  gradient  will  be  much  less.  Since  a  change  of  1°C 
causes  a  change  in  the  refractive  index  of  An/n=10~s.  This 
gives  k=5*10's  km" 1 .  For  /=60°,  and  h= . 5  km,  the  error  in 
the  position  of  the  reflector  will  have  a  component  parallel 
to  the  fault  equal  to  10“4y.  For  y  =  1 0 0  meters,  this  is  1  cm. 
This  error  will  be  reduced  by  about  one  order  of  magnitude 


Figure  8....  Above:  side  view  of  the  path  followed  by 
rays  from  point  P  to  the  camera  at  C.  Below:  top  view  - 
distance  s  split  up  into  its  components  parallel  and 
perpendicular  to  the  fault. 
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if  the  vertical  gradient  of  n  is  approximately  the  same 
during  the  exposure  of  the  second  photograph. 

A  horizontal  temperature  gradient,  if  one  occurs,  will 
be  much  smaller  than  the  vertical  gradient.  A  horizontal 
gradient  of  2°/km  for  example,  at  45°  to  the  line  of  sight 
would  cause  an  error  of  1  mm  in  the  position  of  a  point 
about  1  km  from  the  camera.  It  was  found  by  Vacquier  and 
Whiteman,  (1973)  that  this  error  due  to  horizontal  gradient 
was  too  small  to  affect  their  results  significantly.  Since 
our  lines  of  sight  are  much  shorter  than  theirs,  we  can 
neglect  the  error  due  to  horizontal  refractive  index 
gradient . 

It  is  obvious  that  the  camera  cannot  be  positioned  so 
precisely  that  each  light  would  be  exactly  aligned,  (were 
there  no  movement  along  the  fault  at  all)  and  it  is  also 
unlikely  that  the  average  refractive  index  will  be  the  same 
at  the  time  the  second  photograph  is  taken.  Thus  there  will 
be  a  small  displacement  between  the  lights  as  seen  on  the 
photographs.  In  chapter  4,  it  will  be  shown  that  it  is  not 
only  preferable,  but  necessary  to  add  a  displacement  (not 
necessarily  a  known  displacement)  before  the  comparison  can 
be  made.  Because  of  this  unknown  added  displacement,  only 
non-uniform  strains  can  be  detected.  The  results  of  chapter 
4  also  show  that  high  frequency  noise  does  not  present  a 
serious  problem.  To  get  the  yearly  strain  rate,  a  comparison 
of  two  photographic  records,  separated  in  time  by  one  year, 


will  be  made. 
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A  method  capable  of  detecting  such  small  differences  in 
two  transparencies  will  be  discussed  in  chapter  three, 
several  assumptions  are  made: 

1.  The  same  lens  and  same  camera  were  used  both  times. 

2.  The  lens  was  not  damaged  in  between  the  two  trials. 

3.  The  transparencies  and  their  transforms  were  stored 
carefully . 

4.  The  experiment  was  carried  out  under  conditions  of 
minimum  optical  turbulence  and  under  conditions  of 
similar  temperature  and  pressure. 

5.  Haze  was  removed  by  the  use  of  colour  filters. 

A  photograph  contains  a  vast  amount  of  data,  and 
digitizing  it  would  require  exorbitant  storage  space.  In  our 
case  however,  the  photograph  is  used  to  record  the  precise 
position  of  one  or  two  hundred  small  points,  with  the  rest 
of  the  photograph  being  dark  background.  Optical-digital 
converters  have  reached  a  high  degree  of  sophistication; 
transparencies  can  be  digitized  with  precision  to  500 
lines/mm.  The  experiment  could  therefore  be  equally  well 
carried  out  digitally  or  optically,  (as  will  be  shown  with 
synthetic  data  in  chapter  4) 


2.6  Preliminary  Calculations 

Using  equation  2.13,  and  the  results  of  the  previous 
sections,  it  is  possible  to  calculate  the  width  of  the 
intensity  peaks  on  the  transparency.  It  will  be  shown  in 
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chapter  3  that  the  size  of  the  image  points  ultimately 
limits  the  resolution  of  the  interferometer.  Assume  the  size 
of  the  transparency  is  a  square  decimeter,  and  the  ground 
surface  covers  a  square  kilometer. 

The  following  calculations  are  based  on  parameters  of  a 
commercially  available  terrestrial  camera  (Zeiss  (Jena)  UMK 
10/1318).  The  compound  abberat ion-corrected  lens,  camera 
frame  and  the  film  plate  produce  a  combined  error  no  greater 
than  .003mm.  For  our  purposes,  this  is  not  a  negligible 
error,  however  using  the  same  lens  and  camera  will  eliminate 
most  of  this  error. 

1.  Resolution  of  film  will  be  1000  lines/mm  ,  or  10s 
lines/dm.  This  corresponds  to  an  object  on  the  ground  of 
diameter  1  0 ‘ 5  km ,  or  1cm.  This  is  the  minimum  size  of  an 
object  on  the  ground  which  we  can  resolve. 

2.  Using  concrete  pillars,  it  is  hoped  that  the  camera  will 
not  move  during  exposure,  therefore  the  motion  transfer 
function  ( PM )  equals  one,  even  for  long  exposure  times. 

3.  The  lens  will  have  a  radius  of  2.5  cm  and  a  focal  length 
of  about  50  mm.  Using  equation  2.10,  the  highest  spatial 
frequency  (u)  passed  by  the  lens  is  approximately  1/2/1. 
Taking  the  average  wavelength  of  the  source  to  be  500nm, 
we  get  u= 1 0 6  lines/meter.  This  corresponds  to  an  object 
on  the  ground  with  a  diameter  of  1cm.  Smaller  objects 
will  not  be  passed  by  the  lens.  (The  smallest  object 
passed  by  a  100  mm  lens  of  same  dimensions  would  be  2 
cm)  . 
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4.  For  D=r0=5cm,  we  have  (using  equation  2.20), 

R=.45R+=3500  (cycles /m)  2  (  for  z=1km  ), 
therefore  /  =  .01m  =  1cm. 

These  calculations  pertain  to  long  exposures  (greater 
than  one  second) .  Long  exposures  average  out  very  short 
variations  in  the  index  of  refraction  due  to  turbulence. 
Using  long  exposures  does  not  make  it  unnecessary 
however  to  take  as  many  photographs  as  one  can. 

5.  Since  we  are  taking  an  oblique  photograph,  some  points 
will  be  farther  than  others,  and  thus  some  points  will 
not  be  in  perfect  focus.  The  amount  by  which  a  point  on 
the  ground  is  spread  will  depend  on  the  height  of  the 
camera  and  the  position  of  the  point.  Assuming  the 
points  near  the  central  line  are  in  good  focus,  the 
nearest  and  farthest  points  will  be  spread  by  one  to 
three  cm. 

The  combined  effect  will  be  that  a  10  cm  source  will  be 
imaged  onto  an  area  of  film  larger  than  .01  mm,  but  not 
exceeding  .02  mm.  The  above  errors  affect  the  size  of  the 
image  points,  but  not  their  position.  The  following  errors 
will  alter  the  position  of  the  reflectors. 

1.  The  camera  should  be  positioned  such  that  any  point 
which  has  experienced  no  movement  will  appear  in  the 
same  position  to  within  about  40cm.  Suppose  the  camera 
is  ,3km  above  the  area  to  be  imaged,  and  1km  away  from 
it  horizontally,  then  the  precision  to  which  the  camera 
must  be  positioned  is  7*10's  radians  (vertically),  and 


4*10' 4  radians  (horizontally). 

2.  The  horizontal  gradient  of  n  can  introduce  an  error  of 
about  1  mm  in  the  position  of  the  reflectors.  A  vertical 
gradient  produces  an  error  in  the  displacement  along  the 
fault  which  is  proportional  to  the  distance  of  the  point 
from  the  main  axis.  This  error  will  be  of  the  order  of  1 
mm  for  points  lying  within  100  m  of  the  main  axis. 

3.  Distortions  in  the  size  or  shape  of  the  photographic 
plate  due  to  temperature  or  humidity  can  cause  errors  in 
the  position  of  all  points  in  the  image.  If  we  require 
that  this  error  be  less  than  .1  micron,  (ground 
displacement  error  of  less  than  1  mm)  and  we  are 
measuring  the  distance  between  two  points  which  are 
approximately  .04  mm  apart,  then  the  plate  distortion 
must  be  kept  to  less  than  .25%.  We  can  impose  a  stronger 
condition  on  the  distortion.  If  we  examine  each  pair  of 
dots  by  a  laser  beam  of  1  mm  diameter,  and  we  require 
that  all  points  illuminated  by  the  beam  have  a  position 
error  of  .1  micron,  then  the  distortion  must  be  less 
than  .01%. 

It  is  clearly  demonstrated  by  the  two  photographs  on 
pages  152  and  193  of  Hecht  and  Zajac  (1976),  that  today's 
cameras  are  capable  of  attaining  the  needed  accuracy.  The 
photograph  was  taken  by  a  Hyac  I  panoramic  camera  with  a 
300mm  f/5  lens  (lens  diameter=6cm ) .  It  covers  an  area  of  3 
by  6  km,  and  objects  30  by  60  cm  can  be  clearly  seen  in  the 
blow-up  (second  photograph). 


3.  INTERFEROMETRIC  COMPARISON  OF  TRANSPARENCIES 

Using  a  film  with  resolution  of  1000  lines/mm  enables  us  to 
record  separations  of  10~3mm. 

Each  of  our  two  photographs  will  contain  light  peaks  on 
a  dark  background.  When  superimposed  on  a  single 
transparency  by  a  double  exposure,  the  new  transparency  or 
plate  will  contain  pairs  of  peaks,  separated  by  a  small 
distance  d.  Due  to  the  fact  that  the  camera  cannot  be 
aligned  exactly,  and  that  the  average  index  of  refraction 
was  slightly  different  when  the  first  photograph  was  taken 
and  during  the  exposure  of  the  second  photograph,  d  will  not 
be  entirely  due  to  the  strain  accumulation  on  the  ground 
surface.  It  is  hoped  that  d  will  not  be  more  than  .06  mm, 
which  is  about  three  times  the  width  of  the  peaks,  because 
this  will  limit  the  the  accuracy  of  the  system  (chapter  4). 

Assume  that  the  photographic  plate  has  a  resolution  of 
1000  lines  (or  picture  elements)  per  millimeter,  and  that  it 
has  been  properly  exposed  and  developed.  In  pointwise 
filtering,  we  use  a  laser  beam  as  a  probe  to  examine  a  small 
portion  of  the  photograph  of  diameter  about  1  mm,  which 
corresponds  to  10  meters  of  the  ground  surface  and  thus 
contains  only  one  reflector.  In  whole-field  filtering,  we 
examine  the  whole  photograph  at  once,  but  we  look  at  only 
one  component  of  the  displacement.  Both  methods  will  be 
discussed  in  detail  later. 

From  now  on,  the  units  I  will  use  will  be  picture 
elements.  The  conversion  factor  is  given  by  the  film 


42 


43 


resolution.  In  our  case,  1000  picture  elements  equals  one 
millimeter.  The  diameter  of  the  peaks  thus  becomes  b=20, 
their  separation  will  be  40  to  60,  and  the  laser  beam  covers 
an  area  of  diameter  about  1000. 

Our  method  for  determining  the  displacement  will  be 
similar  to  the  experiment  of  Khetan  and  Chiang  (1976),  and 
Chiang  and  Juang  (1576),  in  which  the  authors  used  laser 
speckles  to  analyse  strain  on  a  plane  surface  under  load. 

The  random  pattern  formed  by  the  speckles  is  a  result  of 
interference  of  coherent  light  bouncing  off  different  parts 
of  an  optically  rough  surface.  The  pattern  is  recorded  on 
high  resolution  film  before  and  after  the  surface  has  been 
deformed  by  an  applied  load,  by  double  exposure  of  the 
photographic  film  or  plate.  The  double  exposure  is  then 
Fourier  transformed,  to  analyse  the  recorded  strain 
information.  The  film  contains  pairs  of  dots,  separated  by  a 
small  distance  d,  one  corresponding  to  the  first  exposure, 
and  the  other  to  the  second.  If  we  consider  a  section  of  the 
film  which  consists  of  only  one  pair  of  dots,  the 
transmittance  of  that  section  approximates  a  pair  of 
pinholes.  The  interference  pattern  produced  by  this  section 
will  therefore  be  a  light  circle  modulated  by  parallel, 
equidistant  fringes.  The  fringe  separation  will  depend  on 
the  separation  and  relative  orientation  of  the  two  pinholes. 

Now,  if  we  examine  the  whole  film,  the  fringes  produced 
by  other  pairs  will  have  different  spacing  and  orientation, 
and  therefore  no  clear  pattern  will  be  visible. 
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This  problem  can  be  overcome  by  two  methods:  Pointwise 
filtering,  or  whole  field  filtering.  In  pointwise  filtering, 
we  examine  an  area  of  the  film  which  is  small  enough  that 
the  displacement  in  that  section  is  uniform.  We  can  do  this 
simply  by  passing  a  laser  beam  (1  to  2  mm  diameter)  through 
the  film,  and  observing  the  diffraction  pattern  on  a  screen. 
In  the  whole  field  method,  we  examine  the  whole  film  at 
once,  and  we  place  a  small  aperture  in  the  frequency  plane. 
We  use  another  lens  to  perform  the  Fourier  transform  of  rhe 
frequency  plane.  The  output  plane  will  then  contain  fringes 
which  are  the  loci  of  points  of  equal  displacement 
component.  By  placing  an  aperture  at  different  positions 
(r,6)  in  the  frequency  plane,  the  fringes  in  the  output 
plane  change  accordingly,  r  changes  the  spacing  of  the 
fringes,  while  9  changes  the  component  of  displacement  we 
are  looking  at.  Thus  the  strain  on  the  test  surface  can  be 
calculated,  however  it  is  useful  to  know  a  priori  if  the 
deformation  was  in  or  out  of  the  plane  of  the  surface. 

Our  problem  is  identical  to  this,  except  that  our 
pattern  of  dots  is  not  completely  random,  and  in  our  case, 
the  pattern  of  dots  is  produced  by  reflectors  positioned  on 
the  ground,  rather  than  a  speckle  pattern  produced  by 
reflection  of  coherent  light  from  a  rough  surface. 

In  the  next  sections,  it  will  be  shown  mathematically 
that  these  fringes  are  a  cosine2  modulation  of  the  light 
intensity  in  the  frequency  plane  of  a  single  exposure.  The 
limitations  which  exisr  on  the  size  of  the  spots  will  also 
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be  given. 


3.1  Recording  a  random  pattern 

In  a  double  exposure  method,  we  expose  each  of  our  two 
photographs  on  a  single  plate,  one  at  a  time. 

Consider  the  complex  amplitude  of  light  h(x1;y,)  in  the 
object  transparency  (the  first  photograph)  at  the  point 
P(x,,y,).  The  light  amplitude  at  P'  due  to  radiation  from  P 
(figure  9)  is  given  by  the  Huygens-Fresnel  principle: 

h ( exp ( i kS ) ) /i\ S  (3.1) 


The  total  amplitude  at  P'  is  therefore: 

•o 

(3.2) 

The  amplitude  g(x3,y3)  at  point  Q  in  the  image  plane  is: 

(see  Appendix) 


[ h ( x , , y , ) /i\ S ] exp( i kS ) dx , ay , 


aO  «0 

PP 


g(x3  ,y3  )  =  jj  [  jj  ( h  (  x  ,  ,  y  ,  ) / i A  S )  exp( i kS )ax , dy , ]A ( x , y ) 


-«0 


•  e 


xp[-ik(x2+y2 )/2f] (exp(ikS'  ) / i A  S ’  )  dxdy 


(3.3) 


where  A  is  the  aperture  function  and  f  is  the  focal  length 
of  the  lens. 


S  2  =p2  +  ( x , -x ) 2  +  ( y , -y ) 2 


Figure  9 


•  •  •  ♦ 


Recording  a  random  pattern 


=  p 2 +  x , 2 +  y , 2  -  2 ( x , x  +y,y)+x2+y2 
=  R 2  -  2 ( x , x+y , y ) +x  2  +y  2 
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S  =  R  "(x,x+y,y)/R  +  (x2+y2)/2R 
=  p  - ( x , x+y , y ) /p  +  ( x  2  +y 2 ) /2p 


(3.4) 


similarly 


S'  =  q  -(xx3+yy3)/q  +  (x2+y2)/2q 


(3.5) 


These  approximations  apply  when  the  plane  P,  and  P3  are  far 
from  the  lens,  and  are  known  as  the  Fraunhofer 
approximations . 

Assume  the  variations  in  S  and  S’  are  small,  then  S=  p 
and  S ’ =q  can  be  considered  constant  (see  Khetan  and  Chiang, 
1976) . 


•exp ( i kS ) exp ( i kS ’ )dx,dy,axdy 


(3.6) 


Now  plug  in  equations  (3.4)  and  (3.5) 


•  exp[ ik (p+q- (xx , +yy , ) /p  -  (xx3  +  yy3)/q 
+(x2+y2)/2p  +  (x2+y2)/2q)]  dxdydx , dy , 

=  -(exp[  ik(p+q)  ]  )  /  /  pq  h(x,  ,y,  )  A(x,y) 

•  exp[(ik/2)  ( x  2  +y  2 ) (  —  1 / f  +  1/p  +1/q  )3 


(3.7) 


•  exp[ - i k ( ( xx , +yy , ) /p  +  ( xx 3 +yy 3 ) /q ) ]  dx , dy , axdy  (3.8) 


out 


1/f  -  1/p  -1/q  =  0 


(3.9) 
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therefore 


g  (  x  3 , y  3 ) =K 


h  (  x  ,  ,  y  ,  )  A  (  x  ,  y  ) 


•  exp[ - i k ( x ( x , /p  +x 3/q ) +y (y , /p  +  y3/q))]  dx , dy , dxdy  (3.10) 
where 


K  =  - (exp[ ik (p+q) ] ) /  \2pq 


Double  exposure  is  made  with  displacement  (as  observed  in 
the  image  plane)  being  d(d*,d^).  Assuming  the  source  shape 
does  not  change  drastically  with  time,  the  total  exposure  at 
the  image  plane  is: 


e(x3 ,y3 )=t[ | g(x3 ,y3 )  |  2  +  | g ( x 3 +d* , v 3 +d; ) | 2  ] 


(3.11) 


where  t  is  the  exposure  time.  Examples  in  chapter  4  show 
that  small  shape  changes  do  not  affect  the  results  very 
much.  We  assume  the  film  was  properly  developed;  then  the 
amplitude  transmission  function  g'  will  be  a  linear  function 
of  exposure  time,  and 


g  ’  (  x  3 , y  3  ) =b  -  ce ( x  3 , y  3 ) 


(3.12) 


The  amplitude  of  light  behind  the  photographic  record  is: 


|g(x3+d, ,y3+ 


g ' (x 3 , y 3 ) =b-ct [ I g ( X 3 , y 3 ) | 2  + 


(3.13) 
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3.2  Whole  field  filtering 

To  obtain  the  displacement  information  from  this 
pattern,  we  perform  an  optical  Fourier  transform  in  the 
usual  way.  We  put  the  transparency  in  the  input  plane,  so 
that  g(x3,y3)  in  the  previous  section  now  becomes  g(x, ,y, ) . 
The  light  intensity  distribution  at  the  transform  plane  of 
the  lens  is  then  the  square  of  the  Fourier  transform  of  the 
amplitude  transmission  function  g'. 


I  ( x  2  ,  y  2  )  = 
•exp[ - i k ( x 


x2 


g’  (x,  ,y ,  ) 

+  y , y  2 )/f ]dx , dy , 


I  2 

i 


(3.14  a) 


Neglecting  the  constant  (or  DC)  term  b,  which  will  give  a 
brignt  spot  at  the  center  of  the  frequency  plane,  we  get 
(away  from  the  center): 


I  (x 

<o 

«f 

JJ 

-  GO 


| g(x , +d* , y , +d ^ ) | 2exp[-ik (x 


exp[ -ik (x  t  x  2+y , y 2 )/f 3dx , dy , 
x  2+y , y 2 )/f  3  ax , dy ,  | 2  (3.15) 


From  the  shift  theorem  of  a  Fourier  transform  (see  Appendix 
for  proof ) ,  we  get : 


P0 

I  (  x  2  ,y  2  )  =c  2 1 2  |  JJ  |  g  ( x  ,  ,  y  ,  )  |  2 exp[  -  i k  ( x  ,  x  2  +y  ,  y  2  ) /f  ] dx  ,  dy  , 

^  ~<C 

+  {  (  | g ( x , , y , ) | 2exp[-ik(x,x2+y1y2)/f]dx1dy1 } 

•  exp[ - i k ( d* x 2 +d^ y 2 ) /f ] | 2  (3.16) 


therefore 
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00 

I  ( X  2  ,  y  2  )  =c  2 1 2  I  (j  I  g  ( X  ,  ,  y  ,  )  I  2  exp[  -  i  k  ( x  ,  x  2  +y  ,  y  2  ) /f  ]  dx  ,  ay  , 

•  {  1  +  exp[-ik(dxx2+a^y2 )/f ] } | 2 

=  I ,  ( x  2  ,  y  2 )  |  1  +  exp[-ik (d„x2+d^y 2 )/f ] |  2  (3.17) 

where 

o 0 

I  i  (  x  2  ,  y  2  )  =c  2 1 2  |  fj  |  g  ( x  ,  ,  y  ,  )  |  2  exp[  -  i  k  ( x  ,  x  2  +y  ,  y  2  ) /f  ]  dx  ,  dy  ,  |  2 

-00 

|  1  +  exp(-iz)  |2  =[  1+  exp(-iz)][  1+  exp(iz)] 

=  1  +  exp(-iz)  +  exp(iz)  +  1 
=2+2  cos  z 

=  2(2  cos 2 ( z/2 )  -  1  +  1  ) 

=  4  cos  2 ( z/2 ) 

therefore 

I ( x  2 , y  2 ) =  4  I , (x2 ,y2 )cos2 [ (k/2f ) (dxx2  +  d^y2)3  (3.18  a) 

where  I,  is  the  intensity  distribution  at  the  transform 
plane  if  only  one  exposure  is  recorded  on  the  film.  Most  of 
the  energy  from  this  term  is  concentrated  near  the  center  of 
the  diffraction  pattern  and  is  referred  to  as  the 
diffraction  halo  of  a  single  exposure  (see  Khetan  &  Chiang, 
1975).  This  halo  is  the  Fourier  transform  of  the  plate  if 
d=0,  or  if  only  one  exposure  has  been  recorded  on  the  plate. 
Its  radius  is  given  by  the  highest  frequencies  in  the 
photograph.  Outside  it,  there  is  no  light  and  therefore  no 
fringes.  It  is  therefore  obvious  that  if  the  photograph 
contains  only  low  frequencies,  we  cannot  observe  fringes  due 
to  small  displacements,  since  they  will  be  outside  the  halo. 
In  other  words,  the  displacement  needs  to  be  larger  than  the 
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width  of  the  intensity  peaks  for  the  whole  field  filtering 
method . 

In  our  case,  the  displacement  is  40,  and  the  peak  width 
is  20.  It  will  be  shown  in  chapter  four  that  there  is  an 
easily  measurable  difference  between  spectra  obtained  when 
d=40,  and  when  d=41.  Suppose  that  the  width  of  the  peak  was 
200  instead  of  20.  Then  we  would  need  to  increase  our 
displacement  to  at  least  200,  and  try  to  measure  a 
difference  between  fringe  patterns  produced  by  d=2Q0  and 
d=20  1  . 

The  resolution  is  thus  limited  by  the  size  of  the 
peaks,  which  is  in  turn  dependent  on  the  atmosphere,  film 
resolution  and  lens  diameter  and  quality. 

When  a  double  exposure  is  made  with  a  displacement 
d(dx  ,dij)  between  exposures,  the  halo  is  modulated  by  cosine 
square  fringes.  These  can  be  seen  if  d  is  uniform  over  the 
whole  field.  If  d  is  non-uniform,  the  corresponding  fringes 
will  have  different  spacings  and  orientations,  and  therefore 
no  pattern  will  be  seen.  If  a  small  aperture  is  placed  in 
the  transform  plane  at  x2,  fringes  will  be  observed  in  the 
image  plane  of  another  lens  depicting  different  values  of 
the  component  of  d  along  x2. 

Dark  fringes  are  obtained  for 

cos[ ( k/2f ) (dxx2+d^y2 ) 3=0 

or  for 


dr  x  2  +dy  y  2  =  (  n  +  1/2  )  A  f 


(3.19) 
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These  fringes  are  the  isothetics,  that  is,  loci  of  points  of 
equal  displacement  component.  For  example,  if  we  are  looking 
at  the  y-component  of  displacement  of  a  disc  under  rotation, 
the  fringes  would  appear  as  straight  lines,  parallel  to  the 
y-axis,  with  spacing  between  fringes  given  by  the  magnitude 
of  the  displacement  component. 


3.3  Pointwise  filtering 

A  selected  point  on  the  photographic  record  is 
illuminated  by  a  narrow  beam  of  coherent  light.  The 
diffracted  light  can  be  observed  by  placing  a  plate  of 
ground  glass  at  some  distance  L  from  the  photographic  plate 
(figure  10).  An  accurate  way  to  measure  these  fringes  would 
be  to  replace  the  ground  glass  by  a  photographic  plate  and 
expose  the  fringe  pattern  directly  on  it.  A 

microdensitometer  could  then  be  used  to  scan  the  fringes  and 
the  obtained  curve  could  then  be  fitted  by  a  least  squares 
method  to  a  curve  of  the  form  A(x)cos2kx,  where  A ( x )  forms 
the  envelope  of  the  fringes.  From  the  geometry  of  the 
figure,  we  have  (see  Khetan  and  Chiang,  1976): 

7  =  { L  2  +  [ x  2 - ( x , -x7  )  ]  2  +  [y  2-  (y  ,  -y",  )  ]  2  } 

=l  [  i  +  (x2+r,-x,  )2/l2  +  (y  2+y7-y ,  ) 2  /l  2  ] 1/2 

Taking  the  first  two  terms  in  the  binomial  expansion, 
(Fresnel  approximation)  we  get: 


F i gur e  1  0 


Arrangement  for  pointwise  filtering 
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/  =  L+  (x22+y22  +  2x2xt  +  2y2y,)/2L 

+  (  (  x  ,  -5T,  )  2  +  ( y  ,  -y”;  )  2  ) /2L  -  (x2x,+  y2y,)/L  (3.20) 


The  intensity  in  the  transform  plane  is: 


x 


y 


( x , ,y, ) exp( i kl ) dx , dy , | 2 


(3.21  ) 


The  first  two  terms  in  equation  (3.20)  will  give  a  phase 
factor,  and  can  be  dropped.  We  can  assume  that  except  near 
the  center  of  the  diffraction  halo, 


( x  2  x , +y  2  y ,  ) /L  >>  ( ( x , -x , ) 2  +  ( y , -y , ) 2 ) /2L  (3.22) 

(This  is  again  the  Fraunhofer  approximation,  therefore  we 
are  again  dealing  with  Fraunhofer  diffraction). 

Inserting  equation  (3.20)  into  equation  ( 3 . 2 1 ) , we  get: 

I  (x  2  ,  y  2  )  =  |  g  '  (x  ,  ,  y  ,  )  exp[  ik  (x  2x  ,  +y  2v  ,  )/L]  dx , dy , | 2  (3.14  b) 

ri»c 

which  is  identical  to  equation  (3.14  a)  except  that  the 
integration  limits  are  given  by  the  beam  size  rather  than 
the  aperture  function.  Then  we  have,  as  before, 

I(x2,y2)=  4  I ,  ( x  2  ,  y  2 )  cos2 [  ( k/2L ) ( dx x 2 +d^ y 2 ) ]  (3.18  b) 

In  this  case,  straight  fringes  can  be  seen  even  for  an 
inhomogeneous  deformation  because  within  the  limits  of  the 
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beam  size,  (  =  1  to  2  mm)  the  displacement  field  can  be 
considered  uniform.  The  fringes  in  this  case  are  uniformly 
spaced  straight  lines,  with  spacing  density  proportional  to 
the  magnitude  of  the  displacement  vector,  and  orientation 
normal  to  the  displacement  vector.  Pointwise  filtering 
yields  more  accurate  results  than  the  whole  field  approach 
for  the  following  reasons: 

1.  The  fringes  are  equidistant  and  their  spacing  can 
therefore  be  averaged. 

2.  The  fringes  are  sharper  and  more  clearly  defined. 

3.  No  lenses  are  required,  thus  there  are  no  errors  due  to 
lens  abberations. 

4.  A  rigid  body  translation  or  rotation  can  be  introduced 
to  give  some  initial  fringes,  and  change  of  fringe 
spacing  can  be  measured,  rather  than  the  fringe  spacing 
itself.  If  no  rigid  body  displacement  is  introduced,  the 
sensitivity  is  the  same  as  for  the  whole  field  case. 

If  we  have  two  peaks,  they  will  be  resolved  (by  Rayleigh's 
criterion)  if  a=b/1.22.  Since  the  intensities  of  the  peaks 
are  added,  we  can't  expect  to  get  any  fringes  if  the  peaks 
are  not  resolved,  therefore 

d ( min ) =  b/1 .22  (3.23) 

Using  the  whole  field  method,  the  smallest  ground 
displacement  we  can  therefore  see  is  20/1.22  =16 
centimeters.  However,  in  the  pointwise  filtering  approach, 
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we  can  add  a  small  displacement  to  the  photographs  before 
they  are  exposed  on  a  single  plate,  to  make  a  about  40.  When 
d=40,  the  Fourier  transform  obtained  has  clearly  defined 
fringes,  with  a  spacing  measurably  different  from  the 
fringes  obtained  when  d=4 1 .  This  means  that  a  difference  in 
ground  displacement  of  1  cm  can  be  detected  using  the 
pointwise  filtering  method,  but  not  the  whole  field  method. 

It  is  possible  to  digitize  the  interesting  sections  of 
the  doubly-exposed  plate,  and  measure  the  separations 
directly.  This  method  would  have  the  advantage  that  if  the 
peaks  were  much  wider  than  the  expected  20  pixels,  (a  pixel 
is  a  picture  element)  and  if  they  were  too  close  to  be 
resolved  by  Rayleigh's  criterion,  they  could  still  be 
deconvolved  (since  we  know  their  approximate  shape),  and 
their  separation  could  be  obtained.  If  the  noise  level  is 
high,  deconvolution  filters  may  become  unstable,  and  even 
this  method  will  not  work.  Assuming  however,  that  it  is 
possible  to  get  peak  widths  of  20  pixels  or  less,  analysing 
the  diffraction  pattern  would  have  the  following  advantages 
over  direct  measurement: 

1.  To  obtain  the  diffraction  pattern  we  need  only  pass  a 
laser  beam  through  the  double  exposure.  The  fringes  can 
then  be  photographed,  and  their  separation  can  be 
measured  with  a  ruler,  in  the  direction  parallel  to  the 
fault  direction.  An  error  in  direction  of  5  degrees  will 
introduce  an  error  of  less  than  1/4  cm  in  the 
displacement  of  the  reflector. 
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To  do  a  direct  measurement  of  peak  separation,  we 
need  to  use  a  microscope  to  enlarge  each  area  and 
digitize  the  enlarged  sections.  This  digitization  needs 
to  be  done  in  two  dimensions.  The  separation  between  the 
peaks  can  be  found  by  cross  correlation  with  a  single 
peak . 

2.  Since  fringes  are  equidistant,  we  can  average  over 

several  fringes  to  get  the  average  separation,  or  do  a 
least  squares  fit  of  the  fringes  to  a  cos2  curve  .  No 
such  averaging  is  possible  when  we  do  a  direct 
measurement . 

If  the  noise  level  is  very  high,  the  direct  measurement 
may  work  better,  since  the  cross-correlation  can  still  pick 
out  the  position  of  the  peaks  even  if  they  are  completely 
buried  in  random  noise.  The  transform  plane  analysis  has 
been  chosen,  because  this  method  makes  it  unnecessary  to 
digitize  in  two  dimensions.  Also,  if  sections  of  the 
original  photograph  are  digitized,  using  a  microscope  and  a 
microdensitometer,  and  the  separations  are  found  using  the 
direct  comparison  method,  we  can  perform  the  digital  Fourier 
transform,  and  analyse  the  digitally-obtained  fringes.  This 
way,  we  will  be  able  to  compare  the  results  obtained  by 
direct  measurement,  and  by  analysing  the  transform  plane. 

In  the  next  chapter,  we  will  assume  that  d=d  and  d  =0, 
thus  enabling  us  to  do  all  computations  in  one  dimension 
without  any  loss  of  generality,  and  at  one  tenth  the  cost 
and  storage  space  requi rememnt .  One  example  is  done  in  two 
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dimensions  to  show  that  small  shape  changes  do  not  affec 
the  results. 


3.4  The  Discrete  Fourier  Transform 

The  calculations  in  chapter  three  were  all  done  by  the 
use  of  the  continuous  Fourier  transforms,  whereas  in  chapter 
four,  the  discrete  Fourier  transform  has  been  calculated  by 
the  computer. 

In  one  dimension,  we  can  use  a  vector  of  10C0  points 
and  perform  the  Fourier  transform  quite  inexpensively.  To 
extend  the  problem  to  two  dimensions,  we  should  first  reduce 
the  number  of  data  points,  since  it  would  be  impractical  to 
use  an  array  of  1000  by  1000  points.  In  order  to  reduce  the 
data,  we  must  first  understand  exactly  how  the  computer 
calculates  the  discrete  transforms.  Normally,  the  Fourier 
series  is  developed  independently  of  the  Fourier  integral. 
Here,  we  are  interested  how  the  discrete  Fourier  transform 
relates  to  the  continuous  one,  therefore  we  will  develop  it 
as  a  special  case  of  the  continuous  Fourier  integral 
(Brigham,  1974  ).  Consider  the  continuous  function  ’n(t),  and 
its  Fourier  transform  H(f),  where  t  can  stand  for  time  or 
space,  and  f  can  be  either  temporal  or  spatial  frequency 
(figure  11).  We  want  a  discrete  Fourier  transform  pair  which 
approximates  the  pair  h(t):H(f).  First,  we  need  to  sample 
the  function  h(t),  by  multiplying  it  by  the  Dirac  comb,  with 
sample  interval  T.  The  transform  of  the  sampled  function 
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Figure  11....  Development  of  the  discrete  Fourier  transform 
pair:  figure  taken  from  Brigham,  1974,  p95 
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h(t)  is  given  by  the  convolution  of  H(f)  and  the  transform 
of  the  Dirac  comb  (see  Appendix  for  proof).  Aliasing  will 
modify  the  transform  pair  (as  seen  in  figure  11(c)),  if  the 
sampling  interval  is  chosen  too  large.  If  h(t)  is 
band-limited,  we  can  avoid  aliasing  error  by  choosing 
T<1/2F,  where  F  is  the  highest  frequency  component  of  H(f). 
So  far,  we  have  only  considered  an  infinite  number  of 

✓A 

samples  of  h(t).  To  be  suitable  for  computation  by  computer, 

/V 

the  function  h(t)  has  to  be  truncated  by  multiplication  with 
a  rectangle  function  x(t),  so  that  N  samples  are  left.  The 
resulting  transform  is  shown  in  figure  11(e).  As  x(t) 
becomes  very  wide,  X(f)  approaches  a  delta  function,  and  the 
ripple  which  has  been  added  to  the  transform  will  be 
attenuated . 

It  is  now  necessary  to  modify  the  frequency  transform 
by  a  frequency  sampling  function  to  make  it  suitable  for 
machine  computation.  Just  as  the  sampling  in  the  time  domain 
resulted  in  a  periodic  function  of  frequency,  sampling  in 
the  frequency  domain  results  in  a  periodic  function  of  time. 
If  the  original  function  is  approximated  by  N  samples,  the 
transform  will  also  be  approximated  by  N  samples  (figure 
11(g)). 

Under  the  following  conditions,  the  discrete  Fourier 
transform  will  be  exactly  equal  to  the  continuous  transform. 

1.  h(t)  must  be  periodic 

2.  h(t)  must  be  band-limited 

3.  Sampling  interval  must  be  smaller  than  1/2F 
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4.  x(t)  must  be  non-zero  over  an  integer  multiple  of 
exactly  one  period  of  h(t) 

The  class  of  functions  we  are  dealing  with  are  of 
finite  duration  in  space,  and  therefore  cannot  be 
band- 1 imi ted ,  but  they  need  not  be  truncated.  Consider  for 
example  a  triangle  function  (figure  12).  In  the  next 
chapter,  the  peaks  will  be  roughly  gaussian  in  shape.  Since 
the  Fourier  transform  of  a  gaussian  is  again  a  gaussian,  a 
triangular  pulse  will  illustrate  these  ideas  more  clearly. 
The  continuous  transform  is  a  Sine2  function.  The  discrete 
transform  may  or  may  not  resemble  the  continuous  one, 
depending  on  the  aliasing  effect,  which  is  determined  by  the 
sampling  interval  T,  and  depending  on  the  spacing  of  the 
points  in  the  frequency  domain,  which  can  be  controlled  by 
the  number  of  zeros  which  we  stick  on  the  ends.  In  figure 
12(a),  N,  the  number  of  samples  is  128,  and  no  zeros  are 
added,  in  figure  12(b),  N  is  16,  with  no  zeros  added,  and  in 
fig.  12(c),  16  samples  are  taken  of  the  triangle,  with  the 

rest  of  the  samples  being  zero,  to  give  a  total  of  128. 
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\  TRANSFORM 


Figure  12....  A  triangle  function  sampled  at  different 
rates,  with  its  discrete  Fourier  transform  on  the  right 


4.  COMPUTATIONS  WITH  SYNTHETIC  DATA 


I  do  not  have  actual  data,  so  in  the  following  computations, 
the  simulated  cross  sections  of  corresponding  sections  of  2 
photographs  were  added  together,  giving  the  intensity 
distribution  of  a  double  exposure  plate.  The  sections  have  a 
diameter  of  1000  picture  elements,  or  1  millimeter,  which  is 
about  the  diameter  of  a  laser  beam  used  to  illuminate  the 
section.  The  double  exposure  is  then  ‘Fourier  transformed,  to 
give  the  intensity  in  the  frequency  plane.  Since  one  picture 
element  on  the  photograph  corresponds  to  1  cm  on  the  ground, 
the  laser  beam  covers  an  area  with  diameter  about  10  meters, 
and  containing  only  one  light.  The  computations  therefore 
simulate  the  pointwise  filtering  approach  described  earlier. 
Maximum  resolution  is  possible  when  the  width  of  the  peak 
(b)  equals  one.  As  the  peak  approaches  a  delta  function,  the 
radius  of  the  halo  approaches  infinity,  and  since  the 
fringes  are  equidistant,  we  can  measure  the  separation 
between  a  large  number  of  fringes,  thus  greatly  reducing  the 
error.  As  the  halo  radius  approaches  infinity,  the  bright 
fringes  approach  a  constant  intensity,  (figure  1) 

Diffraction  patterns  were  obtained  for  d=40,25,10  and 
3,  (figure  2)  and  a  graph  of  (fringe  separation)-1  versus 
the  displacement  d  has  been  made  (figure  3).  It  is  linear, 
as  predicted  by  equation  3.19. 

Next,  we  will  look  at  more  realistic  intensity 
distributions,  that  is,  peaks  of  finite  width.  In  chapter 
two,  preliminary  calculations  showed  that  it  is  possible  to 
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Figure  13....  For  a  peak  of  width  1,  fringes  extend  to 
infinity 


INTENSITY  *10*  INTENSITY  *10’ 

60.00  100.00  160.00  200.00  260.00  cP,*00  60.00  100.00  160.00  200.00  260.00 

_j _ I _ i _ i _ i  . -I - 1 - > - — - 1 - 1 - 1 


65 


INTERFERENCE  D=25 


sit 0.00  560.00  400. 


00 


INTERFERENCE  D=3 


0.00  2^0 .00 

FREQUENCY 


240.00  260 .00  i5o” 


380 . 


Ibo. 


00 


00 


00 


Figure  14....  Displacement  is  inversely  proportional  to  the 
fringe  separation 


66 


Figure  15....  Calibration  graph  showing  inverse 
relationship  between  fringe  separation  and  distance  between 
the  peaks 
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obtain  peaks  of  a  width  of  b=20.  A  displacement  d=0  is 
equivalent  to  a  single  exposure,  and  therefore  shows  only 
the  diffraction  halo  (figure  16).  For  a  delta  function,  a 
displacement  d=2  produces  a  clearly  visible  fringe  pattern. 
In  the  case  of  a  finite-width  peak,  however,  fringes  are  not 
visible,  until  the  displacement  is  larger  than  d(min)  given 
by  equation  3.23.  Since  the  intensities  of  the  two  peaks  are 
added,  we  only  expect  fringes  if  the  total  intensity 
distribution  looks  like  two  peaks,  and  not  one  (figure  17). 
The  diffraction  pattern  produced  by  d=40  is  not  clearly 
different  from  one  obtained  when  d=4 1  or  42,  but  the 
difference  can  be  measured  using  a  ruler  (figure  18). 

The  next  step  is  to  add  random  noise  to  each  photograph 
(figure  19).  The  signal  to  noise  ratio  is  about  10:1  and  the 
fringes  are  still  clearly  visible.  Using  a  digital  band-pass 
filter  subroutine  written  by  Dave  Ganley  (see  Appendix),  the 
high  frequency  noise  was  filtered  out.  In  an  optical  system, 
band-pass  filtering  is  a  very  easy  process,  as  was  shown  in 
chapter  two.  Unfortunately,  the  visibility  of  fringes  was 
not  improved  at  all  (figure  20).  The  same  results  were 
obtained  when  the  signal  to  noise  ratio  was  reduced  to  2:1 
(figure  21)  and  (figure  22). 

The  importance  of  peak  width  has  already  been  discussed 
in  the  previous  chapter.  The  diffraction  halo,  which  forms 
an  envelope  containing  the  fringe  pattern,  decreases  with 
increasing  peak  diameter  (figure  23).  The  peak  width  in  this 
figure  is  about  50. 
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Figure  16....  A  displacement  of  a=0  produces  the 
diffraction  halo,  the  frequency  distribution  due  to 
peak 
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Figure  17....  Fringes  become  visible  only  when  d  >  16,  a 
predicted  by  equation  3.23 
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Figure  18....  Differences  in  diffraction  patterns  for  d=40 
and  d=4 1  can  be  easily  measured  with  a  ruler 
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Figure  19....  Fringes  are  still  clearly  visible  when  random 
noise  is  added  to  the  signal 
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Figure  20....  Amplitude  of  noise  has  been  reduced  by 
band-pass  filtering,  with  no  significant  improvement  in  the 
visibility  of  fringes 
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Figure  22....  At  high  noise  level,  band-pass  filtering 
again  has  no  effect  on  fringe  visibility 
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CROSS  SECTION 


Figure  23....  Radius  of  the  diffraction  hale  decreases  with 
increasing  peak  width.  In  this  figure,  the  peak  is  50  lines 
wide. 
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It  is  to  be  expected  that  the  width  of  the  peak  will 
not  be  the  same  on  two  pictures  which  are  taken  at  different 
times.  If  one  peak  is  wider  than  the  other,  the  size  of  the 
halo  will  be  limited  by  the  size  of  the  wider  peak  (figure 
24).  In  this  figure,  one  peak  had  a  width  of  18  and  the 
second  about  50.  Out  of  the  photographs  taken,  we  hope  to 
have  about  twenty  good  quality  photographs  for  each  year. 
These  will  be  paired  up  such  that  we  are  comparing 
photographs  containing  spots  of  similar  sizes.  Thus 
comparing  two  peaks  with  a  width  ratio  of  3:1  would  occur 
only  in  the  worst  possible  case.  The  error  introduced  by 
this  asymmetry  is  easy  to  calculate  from  figure  24.  The 
fringe  separation  is  inversely  proportional  to  the  distance 
between  the  peaks,  where  the  constant  of  proportionality  is 
32.8  cm2  (calculated  from  figure  13  or  14).  The  mean  fringe 
separation  in  figure  24  is  (.80±.04)cm.  This  gives  a  ground 
displacement  of  41±2  cm;  the  actual  displacement  used  was 
40cm.  The  error  is  large,  not  because  of  the  asymmetry,  but 
because  one  peak  was  three  times  the  maximum  acceptable 
width,  and  because  only  1/2  of  a  diffraction  halo  was 
measured.  In  the  actual  experiment,  we  will  have  about  20 
full  diffraction  patterns,  and  the  error  should  thus 


decrease  by  a  factor 


Clearly,  it  is  undesirable  to  have  two  peaks  in  one 
section,  since  this  will  give  a  superposition  of  two  sets  of 
fringes.  If  the  separation  of  the  lights  is  much  greater 
than  d,  the  superposition  will  resemble  a  beat  pattern, 
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Figure  24  ...  . 
of  the  halo  is 
of  oeak  shown 
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s  50,  other  peak  has 


than  the  other,  the  size 
of  the  wider  peak:  width 
a  width  of  18 
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(figure  25)  and  a  less  enlightening  pattern  is  obtained  if 
the  distance  between  them  is  almost  equal  to  the 
displacement,  or  if  the  displacement  is  different  for  peak. 

If  we  have  a  large  number  of  randomly  spaced  lights  in 
the  area,  then  our  experiment  approaches  that  of  Khetan  and 
Chiang.  The  fringe  pattern  obtained  is  quite  noisy,  and 
resembles  the  pattern  obtained  for  two  peaks,  (figure  26). 
This  is  probably  because  the  pattern  is  not  completely 
random,  but  it  is  obvious  that  an  i nter f erogram  obtained 
with  a  large  number  of  peaks  will  not  have  fringes  as 
clearly  defined  as  that  obtained  using  just  one  peak. 

From  figure  16,  we  know  that  the  maximum  frequency  is 
about  120,  so  the  sampling  interval  used  was  much  finer  than 
necessary.  The  number  of  data  points  can  thus  be  reduced 
from  1024  to  128  before  aliasing  problems  become 
significant.  A  simple  extension  of  figure  18  to  two 
dimensions  has  been  done  (figure  27).  In  this  figure,  the 
width  of  one  peak  is  20  pixels  in  the  y  direction  and  24  in 
the  x  direction.  The  second  peak  is  widest  in  the  direction 
at  45°  to  the  x  axis.  dx=48,  and  d^=40  (|d|=62.48  at  39.8° 
to  the  x  axis).  The  figure  was  drawn  on  an  elect rostat ic 
printer/plotter,  using  16  different  dot  sizes. 
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Figure  25....  A  beat  pattern  is  observed  if  two  lights  are 
present  in  one  section 
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Figure  26....  Interference  pattern  obtained  from  a  random 
pattern  of  delta  functions 
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Figure  27....  Density  plot  of  a  diffraction  pattern 
produced  by  2  gauss ian-shaped  peaks  showing  the  orientation 
and  magnitude  the  displacement.  The  actual  disDlacement  was 
d= ( 48 , 40  )  . 


5.  CONCLUSION 


It  is  only  in  the  last  decade  that  geodetic 
instrumentation  and  measurement  techniques  have  reached  the 
sophistication  necessary  to  detect  yearly  surface 
deformation  near  active  faults.  Stress  accumulation,  which 
is  always  accompanied  by  surface  deformation,  is  of  major 
importance  for  understanding  earthquake  mechanisms. 
Earthquakes  occur  when  the  accumulated  stress  reaches  some 
critical  value,  called  the  failure  stress.  We  proposed  a 
method  of  determining  surface  strain  rates  in  an  area  of 
small  aperture  (about  1  km).  It  was  assumed  that  in  an  area 
10  meters  in  diameter,  the  displacement  will  be  uniform,  and 
the  area  can  therefore  be  represented  by  a  single  point.  A 
row  of  lights  (or  reflectors)  separated  by  10  -  20  meters, 
and  perpendicular  to  the  fault,  can  thus  be  used  to  measure 
the  cross  section  of  yearly  strain  accumulation  across  the 
fault . 

Using  a  rectangular  array  of  reflectors,  rather  than  a 
linear  array  will  increase  the  amount  of  useful  information. 
For  example,  if  the  surface  strain  distribution  is  linear 
across  the  fault,  a  linear  array  of  reflectors  will  still 
look  linear  after  it  has  been  displaced.  One  cannot 
therefore  distinguish  between  a  real  linear  displacement 
distribution,  and  an  apparent  displacement  caused  by 
misalignment  of  the  camera,  or  a  horizontal  temperature 
gradient.  A  rectangular  array,  on  the  other  hand,  will  still 
appear  rectangular  if  the  camera  is  rotated,  however  it  will 
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be  parallelogram-shaped  if  the  displacement  is  real.  Due  to 
the  small  displacements  involved,  I  feel  that  it  is  unlikely 
that  one  could  recognize  the  difference  between  a 
parallelogram  and  a  rotated  rectangle. 

Preliminary  calculations  in  chapter  two  were  based  on 
parameters  set  by  commercially  available  terrestrial  cameras 
or  phototheodolites  (camera  combined  with  a  theodolite). 
Several  models  are  available  from  both  Wild  and  Zeiss 
Canada.  Some  significant  characteristics  of  these  models 
are : 

1.  Picture  formats  vary  from  60  by  80  mm  to  130  by  180  mm. 
on  photographic  plate,  or  cut  and  roll  film. 

2.  Lens  distortion,  camera  frame  and  photographic  plate  or 
film  produce  a  combined  error  no  greater  than  .008  mm. 

3.  Focal  lengths  range  from  45  to  200  mm. 

4.  Focal  length  to  object  distance  should  be  kept  down  to 
1  :  3000  . 

5.  Prices  start  at  $7000  (in  1976) 

It  has  been  shown  that  a  measurable  difference  exists 
in  the  interference  patterns  produced  by  a  reflector  which 
was  displaced  by  40  cm,  and  one  displaced  by  41  cm. 
Non-uniform  displacements  of  the  order  of  1  cm  can  therefore 
be  detected.  To  distinguish  between  real  displacements,  and 
apparent  displacements  caused  by  the  atmosphere,  we  can 
concentrate  on  points  close  to  the  main  axis,  which 
eliminates  error  due  to  vertical  refraction  index  gradient, 
and  we  can  neglect  the  small  effect  of  a  horizontal 
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gradient.  Rapid  atmospheric  fluctuations  are  removed  by  long 
period  averaging,  achieved  by  long  exposure  times  (one 
minute  or  more).  Slow  fluctuations  can  be  minimized  by 
monitoring  the  weather  conditions. 

The  most  important  result  obtained  is  that  by 
introducing  a  small  initial  displacement,  the  sensitivity  of 
the  system  has  been  increased  by  a  factor  of  about  20.  That 
is,  if  the  smallest  detectable  displacement  using  a  direct 
comparison  is  b/1.22,  where  b  is  the  width  of  the  narrowest 
peaks  in  the  photograph,  then  using  the  method  described 
earlier,  the  smallest  detectable  displacement  is  1/20  of 
b/1.22.  At  the  same  time,  the  problems  encountered  in  a 
direct  comparison,  such  as  misalignment,  are  avoided  in  this 
method,  since  the  initial  displacement  can  be  unknown. 


5.1  Suggestion  for  geophysical  experiment 

The  field  set-up  of  the  experiment  is  shown  in  figure 
8.  A  terrestrial  camera  is  used  to  image  an  array  of  lights 
a  distance  h  below.  The  light  sources  could  be  a  set  of 
reflectors,  with  one  powerful  light  source  somewhere  below 
the  camera,  and  well  out  of  sight  from  the  camera's 
position.  Using  reflectors  will  allow  us  to  carry  out  the 
experiment  at  night,  under  less  turbulent  atmospheric 
conditions.  Reflectors  can  be  left  in  the  field  over  long 
periods  of  time,  and  they  will  work  again  when  needed, 
without  having  to  be  located  first.  The  power  needed  for  the 
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single  light  source  will  be  much  greater  than  the  combined 
power  which  would  be  needed  if  all  reflectors  were  replaced 
by  battery-powered  light  sources,  however  due  to  the  large 
number  of  reflectors  proposed  for  this  experiment,  using 
light  sources  would  not  be  economically  feasible. 

The  reflectors  should  be  close  to  the  main  axis,  to 
minimize  the  effect  of  the  vertical  index  of  refraction 
gradient.  They  need  not  form  any  precise  geometrical 
pattern.  A  "random"  pattern  such  as  in  figure  8  is 
sufficient.  The  only  constraint  is  that  the  reflectors  lie 
in  a  band  approximately  200  meters  wide,  centered  on  the 
main  axis.  The  ultimate  goal  of  the  experiment  is  to  overlap 
two  photographs  of  these  lights  onto  a  single  transparency. 
This  will  produce  an  array  of  light-pairs  on  a  dark 
background.  The  separation  and  relative  orientation  of  the 
two  lights  in  each  pair  can  then  be  determined  by  examining 
the  interference  pattern  they  produce  when  a  laser  beam  is 
passed  through  them. 

The  experiment  must  be  broken  up  into  sections,  and 
each  part  must  be  carefully  tested  before  going  on  to  the 
next.  In  the  first  case,  the  width  of  each  light  on  the  film 
must  be  kept  to  less  than  .02  mm.  Theoretically,  this  can  be 
easily  achieved  using  a  film  with  resolution  of  1000 
lines/mm,  and  a  f/1.2  lens.  The  light  source  will  produce  a 
bright  spot  on  the  transparency.  The  width  of  the  spot  can 
be  easily  obtained  by  inspecting  the  Fraunhofer  diffraction 
produced  when  we  pass  a  laser  beam  through  the  spot.  The 
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film  must  be  properly  developed  to  get  peaks  of  this  width. 

This  test  should  be  repeated  now  with  a  large  distance 
between  the  light  source  and  camera,  to  simulate  the 
atmospheric  conditions.  The  required  peak  width  should  again 
be  less  than  .02  mm.  Here,  it  is  important  to  have  a  lens 
with  diameter  of  5  cm  or  more. 

Once  the  double  exposure  has  been  obtained,  and  the 
separations  between  the  peaks  in  every  pair  have  been 
measured,  the  separations  can  be  graphed  as  a  function  of 
the  peaks'  perpendicular  distance  from  the  fault.  Before 
graphing  the  results,  we  must  remove  the  displacement  due  to 
camera  alignment.  To  do  this,  we  will  assume  that  the 
closest  reflector  has  not  moved  with  respect  to  the  camera, 
and  that  the  camera  is  in  exactly  the  same  position,  but  may 
have  been  rotated  slightly. 

Using  the  model  of  Turcotte  and  Spence  (1974,  equation 
6)  one  could  do  a  least  squares  fit  to  the  data,  to  find  the 
rate  of  strain  accumulation,  and  the  depth  below  which  the 
plates  slip  freely. 

The  errors  in  the  position  of  the  reflector  images  are 
due  to  the  atmosphere,  camera  alignment,  and  plate  stability 
limitations.  The  errors  due  to  the  atmosphere  are  most 
troublesome,  because  of  their  variability.  As  was  found  by 
Vacquier  and  Whiteman,  a  large  percentage  of  the  photographs 
are  unusable  because  of  unfavorable  atmospheric  conditions 
during  exposure. 
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Photographic  plates  should  be  stable  to  within  .01%.  A 
10  cm  plate  should  therefore  be  stable  to  10  microns. 

Testing  plate  stability  can  be  done  in  the  lab,  as  the  first 
step  in  determining  the  practical  feasibility  of  this 
method . 

The  camera  cannot  be  aligned  precisely  as  it  was  one  or 
two  years  before.  This  unfortunately  limits  our  method  to 
detection  of  non-linear  displacement  fields. 

The  ideas  behind  this  experiment  are  very  simple,  but 
to  actually  cary  out  the  whole  experiment  will  require  a 
great  deal  of  care  and  patience.  The  method  seems  feasible 
theoretically,  and  this  work  should  be  taken  into  its  next 
stage,  which  is  to  determine  the  practical  feasibility. 
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APPENDIX  A:  Notation 


f=focal  length 

A=wavelength  of  source  (=2nc/w) 

P,=object  plane 
P2=transform  plane 
P3=image  plane 
( x , , y , ) =point  in  P, 

(x2 ,  y 2 )=point  in  P2 
( x 3 , y 3 ) =point  in  P3,  inverted 

( u , v ) =spat ial  frequencies  in  x  and  y  directions  respectively 
x  =  ( x  ,  y )  ,  y  =  ( u ,  v ) 


i (x ,  y )= input 
o ( x  ,  y ) =output 

p(x,y)=point  spread  function  (impulse  response) 

Fourier  transforms  of  functions  are  denoted  by  capital 
letter  s 

? ( u , v ) = opt i cal  transfer  function  (F.T.  of  point  spread) 
s ( x , y ) =observed  image 
n ( x , v ) =no  i  se 

A ( x , v ) =apert ure  function 

I,0,P,S,N  are  Fourier  transforms  of  i,o,p,s,n  respectively 
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APPENDIX  B:  Proofs  of  Fourier  Transform  Theorems 

Fourier  Transform  properties  of  lenses 

Consider  the  input  amplitude  distribution  g, (x, ,y, )  in 
front  of  a  spherical  lens,  with  the  observation  plane  a 
distance  f  behind  the  lens(see  figure  5).  From  the 
Huygens-Fresnel  principle,  we  have: 

g(x,y)  = [exp( ikd)/iAd] 

•  jjg  ,  exp[  i  k/2d  (  (  x-x  ,  )  2  +  ( y-y  ,  )  2  )  ]  dx  ,  dy  ,  A1 

where  we  used  the  approximations 

c  o  s  ( n  ,  r  ,  )  =  1 

r,=  d[  1  +  (x-x,)2/2d2  +  (y-y,)2/2d2] 

and  we  ignore  the  finite  size  of  the  aperture.  If  the 
constant  phase  factor  is  dropped,  then  we  have 

g(x ,y )=( 1/iAd)  jj  g,(x,,y,)  p( x-x , , y-y , ) dx , dy , 

=exp[ik/2d  ( x 2 +y 2 ) ] / i\ a 

•  jj  g  ,  exp[  i  k  (  (  x  ,  2  +y  ,  2  ) /2d  -  (  x  ,  x+y  ,  y  )  /d )  ]  dx  ,  dy  ,  A2 

similarly, 

g2  (xa  ,y2  )  =exp[  (  i  k/2f  )  (  x  2  2  +y  2  2  )  ] /U  f 

•JJg'  exp[  i  k  (  (  x  2  +y  2  ) /2  f  -  {  x  2  x+y  2  y  ) /f  )  ]  dxdy  A3 


substitute  in  for  g'  from  equation  2.5 
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A4 

AS 


A6 


A7 

Obviously,  if  d=f,  the  first  exponential  drops  out  and  we 
are  left  with  an  exact  Fourier  transform  relation. 

Shift  theorem 

F{g(x-a,y-b)}  =  g  (x-a  ,  y-b)  exp[-2ni(ux  +  vy)  ]dxdy 
=  jfg  (  x  ’  ,  y  '  )exp[-2ni(u(x'+a)+v(y'+b))]dx’dy' 
g ( x ’  , y ’  ) exp[ -2n i ( ux  f +vy ' ) ]dx ' av ' 

•exp[-2ni (ua+vb) ] 

=G(u,v)  exp[-2ni (ua+vb) ]  A8 

Convolution  theorem 

I  f  G ( u , v )  =  F{  g(x,y )  } 
and  H ( u , v )  =  F{  h(x,y)  } 
then 

F  {  g  (  x  '  ,  y  '  )  h  (  x - x  ’  ,  y - y  ’  )  dx  '  dy  '  } 


g  2 ( X  2 , y  2 )  =exp[ ( i k/2  f )  (x2  2+y2  2 ) ] / i A f 

x,y)exp[ ( — i k/f ) (x2x+y2y) ]dxdy 
=  (exp[ ( i k/2  f ) ( x  2 2+y2  2 ) ]/i\ f )  G(u,v) 
where  ( u  ,  v  )  =  (  x  2  /\  f  ,  y  2  /\  f  ) 

G ( u , v ) =G , ( u , v ) P ( u , v )  (from  equation  A1 

P(u,v)=exp( ikd)exp[-inAd(u2+v2 ) ] 

=  exp(ikd)exp[ (-ikd/2f  2 ) (x2  2+y 2  2 ) ] 


G ,  = 

.heref ore 


oo 


g ,  ( x  ,  , y  ,  )  exp[ ( -ik/f ) ( x , x  2  +y , y  2 ) ] dx , dy , 


g j =exp[ ( i k/2f ) ( x 2 2 +y 2 2 ) ( 1  -  d/f ) ] 

•  tv  g  ,  ( x  ,  ,  y  ,  )  exp[  ( -ik/f  )  (x  ,  x  2  +  y,y2)dx,dy, 
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g(x',v')  F{h(x-x ' , y-y ' ) }  dx ’ dy ' 
g ( x ' , v ' ) exp[ -2n i ( ux ' +vy ' ) ]  dx ' dy '  H(u,v) 
(u ,  v )  H ( u r  V ) 


A9 


2 

3 

4 

5 

6 

7 

8 

Q 

1  o 

1  1 

1  2 

1  3 

1  4 

1  5 

1  6 

1  7 

1  8 

20 

2  1 

22 

23 

2  4 

25 

26 

27 

2  8 

29 

30 

3  1 

32 

33 

3  4 

34 

35 

36 

37 

38 

3  9 

40 

4  1 

42 

4  3 

4  4 

4  5 

46 

4  7 

4  8 

4  9 

50 

5  1 

5  2 

53 

54 

55 

56 

57 

5  8 

5  9 

6  O 

6  1 

6  2 

6  3 

6  4 

6  5 

6  6 

6  7 

S  8 

6  9 

70 

7  1 

72 

87 

88 

8  9 

90 

9  1 

9  2 

93 

9  4 

9  5 

9  6 

9  7 

9  8 

9  9 

1  OO 

1  O  1 

102 

1  03 

1  04 

1  05 

1  06 

1  07 

1  08 

1  09 

1  1  O 

1  1  1 

1  1  2 

1  1  3 

1  1  4 

1  1  5 

1  1  6 

1  1  7 

1  1  8 

1  1  9 

1  20 


98 


APPENDIX  C:  Programs 


DIFFRACTION 


BY  GEORGE  KOVflR  JANUARY  1982 

THIS  PROGRAM  USES  THE  IMSl  FAST  FOURIER  TRANSFORM 
TO  OBTAIN  THE  DIFFRACTION  PATTERN  OBSERVED  WHEN  A 
LASER  BEAM  PASSES  THROUGH  A  PAIR  OF  PIN  HOLES 
DIMENSION  IWK  (  1001  ,  FTCSOl  1024)  ,WK(  100)  ,  FTCS01  <  1024 

‘  .  D  (  8  ) 

COMPLEX  C ( 1 024 ) 

N:  1  02  4 
JN  =  4  1 

DISPLACEMENT  IN  FRAME  UNITS  IS  EOUAL  TO  JN-1 
R  E  A  D  I  8  ,  8  8  )  I  FTCSOl  J  )  .  J  =  1  ,N)  ,  (FTCSOl  <J)  ,J=1  ,N> 

FORMAT  (  1  X  ,  1  6  F  4  0  ) 


90  FORMAT ( 1 X , 8F9 . 3 ) 


CALL 

PLOTS 

CALL 

P  L  OTTR  (  FTCSO  . 

N  ) 

CALL 

SYMBOL (1,5.1 

,  ■  3 

,  'CROSS  SECTION 

CALL 

BNOPAS  (  30 . 0  ,  1 

OO 

,  1  .  , 

,0,0) 

CALL 

F  I  LTER  (  FTCSO , 

N  ,  D 

,  G  ,  ‘ 

1  ) 

CALL 

F I  L  TER  (  FTCSO 1 

.  N  , 

D  ,  G  , 

,  1  ) 

CALL 

P  L  0  T  (  -  2  .  ,6  ,  - 

3  > 

CALL 

PL  OTTR  (  FTCSO  , 

N  ) 

CALL 

SYMBOL  <1,51 

,  ■  3 

.  '  F  1 

1  L  TERED  CROSS 

00  5 

JsJN, N 

FTCSO (  J  )  =  F  T  C  S  0 (  J  )  - 

FTCSO 1 i 

(  J  -  JN+  1  ) 

DO  6 

J  =  1  ,  N 

FTCSO  t  J  >  =  ABS  (  FTCSO 

<  J  ) 

> 

Ci  Ji  = 

FTCSO  (  J  ) 

WRITE 

(6,90)  (FTCSOl 

J  >  . 

J  =  1  , 

,  N  ) 

CALL 

FFTCCI  C,  N,  IWK 

,  WK 

> 

00  3 

J  =  1  ,  N 

FTCSO 

(  J  )  =  CA8S  (  C  (  J  ) 

t  *  * 

2 

WRITE!  6  .  90)  (FTCSOl 

J  )  . 

J  =  1  , 

,  N  » 

CALL 

PLOT  (  -  2  .  ,6.  ,  - 

3  ) 

N  =  N  /  5 

CALL 

PLOTTR ( FTCSO, 

N  ) 

CALL 

SYMBOL  <2,51 

,  •  3 

.'INTERFERENCE  D 

CALL  PLOT  I  O  .  ,  O  .  ,  999  ) 

RETURN 

END 

SUBROUTINE  PLOTTR(Y.N) 

BY  G.  KOVAR,  JANUARY  1982 
DIMENSION  Y  (  1 026  )  .  X  (  1 026  ) 

THE  X  COMPONENT  IS  INCREASED  BY  A  CONSTANT  INCREMENT 
THE  Y  COMPONENT  COMES  FROM  THE  CALLING  PROGRAM 
STRAIGHT  LINES  ARE  ORAWN  BETWEEN  ADJACENT  POINTS 
N  IS  THE  NO  OF  DATA  POINTS 


FORMAT (IX, 
00  10  1=1, 
X  (  I  )  =  I 


8  E  9 
N 


3  > 


CALL  • 

P  L  OT  <  2  .  , 

1  .  .  - 

3  ) 

CALL 

SCAL  E (  Y  , 

5  .  ,  N 

.  1  > 

CALL 

SCAL  E  (  X  . 

1  O  , 

N  ,  1  ) 

CALL 

AX  I S  (  0  . 

0  .  ,  ' 

I  NTENS  ITY  ' 

,9.5 

.  90 

Y  ( 

IF  (  N 

.  LT  .  1 OOO  )  GO 

TO  75 

CALL 

AX  I S  (  0  .  , 

0  .  ,  ' 

DISTANCE 

,  -  9  , 

10.  ,  0  . 

,  x 

GO  TO 

76 

75 

CALL 

AX  I S  (  0  , 

0  .  .  ' 

FREOUENCY ' 

,  -  9  , 

10  .0 

,  X 

7  6 

CALL 

F  L  I NE  (  X  , 

Y  ,  N  , 

1  ,0,0' 

N  2  =  N  +  2 

WR  I  T  E  (  6  ,  9  0  )  (  Y 

(  I  )  , 

1=1  ,  N2  ) 

•  1  )  ,  Y  (  N  +  2  )  ) 


X(N*1  )  ,  X  (  N  +  2  )  > 


X  (  N  +  1  )  ,  X  (  N  +  2  )  ) 


RETURN 

ENO 


SUBROUTINE  8N0PAS(F1  ,F2,DELT,D,G) 

BY  OAVE  GANLEY  ON  MARCH  5,  1977 

THIS  SECTION  CALCULATES  THE  FILTER  AND  MUST  8E  CALLED  BEFORE 
FILTER  IS  CALLED 


FI  =  LOW  FREQUENCY  CUTOFF  (6  DB  DOWN) 

F  2  =  HIGH  FREOUENCY  CUTOFF  (6  OB  DOWN  I 

OELT  =  SAMPLE  INTERVAL  IN  MILLISECONOS 

D  =  WILL  CONTAIN  8  2  DOMAIN  COEFICIENTS  OF  RECURSIVE  FILTE 
G  =  WILL  CONTAIN  THE  GAIN  OF  THE  FILTER 

WRITE  (6,1)  FI .F2.DELT 

FORMAT  ('1  BANOPASS  FILTER  DESIGN  FOR  A  BANO  FROM  ',F8.3,'  TO 
3.'  HERT2  ',//'  SAMPLE  INTERVAL  IS  ',F5.2,'  MILLISECONDS  ') 
DT=0ELT/1 OOO . O 
TDT  =  2  .  0/OT 
F  0  T  =  4  O/DT 
I  SW=  1 


P(  1  >  =  CMP  L  X  (  -  .3826834, 

P(2)=CMPLX(  -  .  3826834  , 

P  (  3  >  =  CMP  L  X  (  -  .  9  238795  , 

P(4)=CMPLX(  -  .  9  2  3  8  7  9  5  . 

W 1  =  TWOP  I  *  F  1 

W2=  TWOP I  *  F2 

W1  =  TUT  *  TAN(W1/TDT) 

W2  =  TOT* TAN (W2/TDT  ) 
HWlDs IW2-W1 ) / 2 . O 
WWsWI *W2 
00  19  1=1,4 

2  1 =  P  (  I  >  *  HWI D 
22=21*21 -WW 
2  2  =  CS  OR  T (  22  > 

S  (  I  )  =  2  1  +22 
S< 1+41=21-22 
WRITE  (6.2)  S 
FORMAT  ( ' - S  PLANE 


9238795  ) 

■  .  9238795  > 
3826  834  ) 

•  3826834  > 


POLES  ARE  AT 


,  8  (  /  ' 


1  8  1 

1  8  2 

1  8  3 

1  8  a 

1  85 

1  8  6 

1  8  7 

1  8  8 

1  8  9 

1  90 

1  9  1 

1  9  2 

1  93 

1  9  4 

1  9  5 

1  9  6 

1  97 

1  9  8 

1  9  9 

200 

20  1 

202 

203 

204 

205 

206 

207 

208 

209 

2  1  0 

2  1  1 

2  1  2 

2  1  3 

2  1  4 

2  1  5 

2  1  6 

2  1  7 

2  1  8 

2  1  9 

2  20 

22  1 

222 

223 

22  4 

2  2  S 

226 

227 

226 

229 

230 

1  2  1 

1  22 

1  23 

»  24 

1  25 

1  26 

1  27 

1  2  8 

1  29 

1  30 

1  3  1 

1  32 

1  39 

1  40 

1  4  1 

1  4  2 

1  43 

1  4  4 

1  45 

1  46 

1  47 

1  48 

1  4  9 

1  SO 

1  5  1 

1  5  2 

1  S3 

1  54 

1  5  5 

1  56 

1  57 

1  5  8 

1  5  9 

I  60 

1  6  1 

1  6  2 

1  6  3 

1  6  4 

1  6  5 

1  6  6 

1  6  7 

1  6  8 

1  6  9 

1  70 

t  7  1 

1  72 

1  73 

1  7  4 

1  75 

1  76 

1  77 

1  78 

1  79 

1  80 


99 


C 


c 

c 

c 

c 

c 

c 

c 


c 


3  6  Ms  3 

M  1  =  2 
M  2  =  1 

37  XCIMI  =  XM-XM2-0I  1  >*XClM1  '  *  0  (  2  l  « X C  ( M2  ) 

XD ( M I sXC I M I  *  XC I  M2  I  -0<  3  i *  XO ( M 1  )  -  0  (  4  ) *  X  D ( M  2  > 

XE  ( M  I : XO  l M )  -  X  D  <  M  2  >  -0<  5  I  *  XE I M  I  )  -  D  <  6  )  *  XE  <  M2  > 

39  X(  I)rXE(M)-XEIM2)-OI7)*X(  I-n-D(8)*Xl  1-2) 

FILTER  IN  REVERSE  DIRECTION 
X  M  2  =  X  i  N  ) 

XMUXIN*  I  I 
XM: X  (  N  -  2  > 

XC  (  1  >  s  X M 2 

X  C  (  2  )  *  XM  »  - D (  1  >  *  X  C  t  1  ) 

X  C  (  3  I  s  XM • XM2 • D I  1  ) *  X C  I  2  I  • 0 (  2  I  *  X C  I  1  ) 

X0(  1  )  =  XC (  1  ) 

XD(2)  =  XCI2)*0(3I«XD(  1  ) 

X0<3)sXC(3)-XCl  1  )-D(3)*XD<2>-0<4)«X0(  1  > 

XE(  1  ):X0(  1  ) 

XE ( 2  )  s  XO (  2)  -  0(5  >  *  X  E  (  1  ) 

XE(3)  =  X0(3)*X0I  1  >*0<5)*XE(2)-0(6  I  *  XE  (  1  ■ 

X  (  N  )  *  XE (  1  ) 

X  <  N  *  1  )cXEI2)  -  01  7I*XI  1  ) 

X(N*2)sXE'3)*XE(  1  l-DI7>*X(2)*0i  8>*X(  1  > 

DO  49  1  =  4  ,  N 

X  M  2  = XM 1 
XM  1  s  XM 
JsN-I+1 

XM  s X  (  J  f 

K*  I  -  (  (  I  -  1  >  /3  >  «3 
CO  TO  (  44  ,  45  ,  46  )  ,  K 

44  Ms  1 
M  1  =  3 
M2  =  2 

GO  TO  47 

45  Ms  2 
M  1  s  1 
M  2  =  3 

GO  TO  47 

4  6  Ms  3 

M  1  s  2 
M2  =  1 

47  XC  ( M  )  s  XM*  XM2  •  0  I  1  >*XC<M1  )  -  0 (  2  )  *  X C  ( M2  I 

X0IM|rXCIM)-XC(M2)*0(3)*X0(M1 ) -0(4I*X0(M2) 

XElMlsX0(M)*X0'M2)  *0(5  '  *  XE  (M  I  )  -  0  (  6  I  *  X  E ( M2  ) 

49  X<  J  I  s  XE I M )  • XE { M2  )  “0(7  )  *  X  (  J+1  )  *0  I  8  I  « X  <  J+2  ) 

IF  <  I  G  .  NE  .  1  )  RETURN 
DO  59  I  s  1  ,  N 
59  X  <  I  )  s  x  (  I  )  /  G 
RETURN 
ENO 

G=  5 / HWI D 
G  =  G  *  G 
Gs  G*  C 

DO  29  Is  i  ,  7,2 
8=  -  2  0*REAL  (  S  (  I  )  > 

2 1 sS <  I  ) *S (  I  ♦  1  ) 

CsRE AL  (  2  1  ) 

AsTDT+e+C/TDT 
Gs  G*  A 

D  <  I  )*(C*DT-FOT)/A 
29  0( I+t )s(A-2.0*8)/A 

Gs  G*  G 
RETURN 

ENTRY  FILTERIX,N,D,C, IGl 

X  =  OATA  VECTOR  OF  LENGTH  N  CONTAINING  DATA  TO  BE  FILTERED 
0  =  FILTER  COEFFICIENTS  CALCULATED  3Y  BNDPAS 
G  =  FILTER  GAIN 

IG  s  1  MEANS  TO  REMOVE  THE  FILTER  GAIN  SO  THAT  THE  GAIN  IS  UNITY 

IF  I  ISW.EO.  1  I  GO  TO  3  1 
WRITE  (6,6) 

6  FORMAT  (  '  1  BNOPAS  MUST  BE  CALLED  8EF0RE  FILTER  - 
CALL  EXIT 

APPLY  FILTER  IN  FORWARD  DIRECTION 
3  1  X  M  2  s  X  M  ) 

XMUXI  2) 

XM: X ( 3 ) 

XC(  1  >  =  XM2 


XC  (  2 

) 

s  XM 

1 

-  D  (  1  >  » 

X  C  (  1  » 

X  C  (  3 

) 

s  XM 

- 

X  M  2  -  0  ( 

!  I*XC 

t  2 

>  -  0  (  2  )  * 

XC 

'  1  * 

XO  t  1 

) 

=  X  C 

( 

1  ) 

XO  (  2 

) 

sXC 

< 

2  »  -  D  (  3 

)  «  XO  t 

]  i 

X  D  (  3 

) 

sXC 

( 

3  )  -  X  C  i 

1  )  -  D  ( 

3  ) 

*  X  D  <  2  )  - 

D  i 

4  )  * 

X  D 

X  E  (  1 

) 

s  X  D 

i 

1  ) 

X  E  (  2 

> 

s  X  D 

( 

2)  *0(5 

)  *  XE  ( 

1  ) 

XE  (  3 

) 

s  XD 

( 

3  )  -  X  D  < 

1  )  -  D  ( 

5  I 

*  X  E  <  2  >  - 

0  ( 

6  >  * 

XE 

X  <  1  > 

X 

XE  ( 

1 

) 

X  <  2  ) 

z 

XE  ( 

2 

)  -  0  (  7  ) 

*  X  <  1  ) 

X  (  3  ) 

S  X  E  ( 

3 

>  -  XE  (  1 

>•0(7 

'  * 

X  (  2  )  -  0  ( 

8  > 

•  X  ( 

1  i 

00  39 

I  = 

4 

,  N 

XM 2  *  KM  1 
XM 1  s  XM 
XM:  X  <  I  ) 

K  =  I  *  (  (  1-1  )  /  3  >  *  3 
GO  TO  (34,35,36)  ,  K 

34  Ms  1 

M  1  s  3 
M2=  2 

GO  TO  37 

35  Ms  2 

M  1  s  1 

M2  =  3 
GO  TO 


37 


